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ABSTRACT:  The  NOL  Hypervelocity  Wind  Tunnel  will  provide  a  high 

Reynolds  number  turbulent  flow  simulation  in  the  Mach  number  range  10 
to  20.  This  facility,  much  needed  for  large-scale  testing  of  hyper¬ 
sonic  vehicles,  is  under  construction  and  will  be  operational  late  in 
1972.  Supply  pressures  up  to  40,000  psi  will  be  maintained  constant 
for  1  to  4  seconds  during  which  stable,  condensation-free  flow  con¬ 
ditions  will  prevail.  Very  high  Reynolds  numbers,  ranging  from  6.5  x 
10®  at  Mach  number  20  to  46  x  10®  at  Mach  number  10  (based  on  the 
nozzle  exit  diameter) ,  are  obtained  by  operating  with  nitrogen  at 
temperatures  just  sufficient  to  avoid  test  section  condensation.  This 
report  presents  the  detailed  procedure  used  to  design  three  nozzles 
for  this  facility.  The  nozzles  are  forty  feet  long  with  five  foot 
exit  diameters  and  are  designed  to  operate  at  Mach  numbers  10,  15,  and 
20  with  supply  pressure  of  430,  2365,  and  3110  atmospheres,  respectively. 
Pressure  and  temperature  effects  on  the  thermodynamic  properties  of 
nitrogen  are  important  at  these  elevated  supply  conditions  and  are 
taken  into  account  in  both  the  inviscid  core  and  the  boundary  layer 
calculations.  The  supersonic  portion  of  the  inviscid  core  is  cal¬ 
culated  by  the  method  of  characteristics  and  the  subsonic  portion  is 
approximated  as  a  one-dimensional  flow.  The  boundary  layer  calcula¬ 
tion  procedure  is  based  on  an  integral  moment- of- momentum  equation 
and  is  valid  for  thick  as  well  as  thin  boundary  layers.  The  design 
conditions  and  the  coordinates  of  both  the  nozzle  wall  and  the  inviscid 
core  contours  are  given. 
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This  report  is  one  of  a  series  of  technical  reports  describing  the 
design  and  performance  of  the  NOL  Hypervelocity  Wind  Tunnel.  A 
previous  report  describes  the  facility  and  its  mode  of  operation  in 
very  general  terms  and  summarizes  the  more  significant  aspects  of  the 
aerodynamic  design  of  the  facility.  The  present  report  describes  more 
completely  the  detailed  procedure  used  to  design  the  nozzles  for  this 
facility.  Future  reports  in  the  series  will  include  more  detailed 
descriptions  of  individual  design  features  or  calculation  methods  used 
in  the  development  of  the  facility.  Tunnel  performance  and  test 
results  will  be  reported  when  available. 
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INTRODUCTION 

A  new  hypervelocity  wind  tunnel  is  currently  being  installed  at 
the  U.S.  Naval  Ordnance  Laboratory  (NOL) .  The  facility,  which  operates 
for  a  minimum  of  one  second  in  a  blowdown  mode,  has  separate  legs  for 
Mach  number  10,  15,  and  20  operation.  Each  leg  is  fitted  with  a  storage 
heater  and  a  forty-foot  long,  five-foot  exit  diameter  contoured  axisym- 
metric  nozzle,  but  shares  a  common  test  section.  Very  high  Reynolds 
numbers,  ranging  from  6.5x10  to  46x10  based  on  nozzle  exit  diameter, 
are  obtained  by  operating  with  nitrogen  at  pressures  up  to  3110  atm  and 
temperatures  just  sufficient  to  avoid  test  section  condensation. 

Various  aspects  of  the  mechanical  and  aerodynamic  design  of  this 
facility  were  summarized  in  reference  (1) .  The  aerodynamic  design  was 
presented  in  somewhat  more  detail  in  reference  C21 .  This  report  will 
present  more  completely  the  detailed  procedure  used  to  design  the  three 
nozzles  for  this  facility. 


PROPERTIES  OF  HIGH-DENSITY  NITROGEN 
THERMODYNAMIC  PROPERTIES 

The  expansion  of  a  gas  in  a  nozzle  is  an  isentropic  process, 
provided  the  gas  remains  in  local  thermodynamic  equilibrium  during  the 
expansion.  If  the  expansion  is  not  too  fast,  the  small  departures 
from  local  thermodynamic  equilibrium  can  be  neglected  and  the  flow 
considered  to  be  isentropic.  The  gas  properties  required  for  the 
nozzle  design  calculations  are  then  the  equilibrium  thermodynamic 
properties  of  the  gas  and  as  such  can  be  tabulated  once  the  desired 
isentrope  is  identified.  In  practice,  the  isentrope  can  be  determined 
by  choosing  either  the  nozzle  supply  conditions  or  the  test  section 
conditions  desired.  Once  the  required  gas  properties  are  tabulated, 
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only  the  tabulations  need  be  introduced  into  the  nozzle  contour  cal¬ 
culations.  This  in  effect  uncouples  the  properties  calculation  from 
the  contour  calculation  and  thereby  simplifies  the  nozzle  design 
problem  considerably.  For  the  present  nozzles,  the  flow  conditions 
and  nozzle  sizes  discussed  in  later  sections  indicate  that  the  nitrogen 
will  be  very  close  to  local  thermodynamic  equilibrium  throughout  the 
expansion  process.  The  small  amount  of  energy  frozen  in  the  vibrational 
mode  will  have  little  effect  on  the  overall  properties  and  can  be 
neglected.  Thus,  for  the  three  nozzles  discussed  in  this  report,  the 
nitrogen  test  gas  is  assumed  to  remain  in  local  thermodynamic  equilibrium. 

In  a  later  section,  it  will  be  shown  that  in  order  to  achieve  the 
desired  range  of  free-stream  Reynolds  numbers  in  the  tunnel  test- 
section  at  reasonable  mass  flow  rates,  it  is  necessary  to  use  nozzle 
supply  pressures  of  several  thousand  atmospheres  for  operation  at  Mach 
numbers  of  15  and  20.  The  corresponding  densities  are  sufficiently 
high  that  compressibility  corrections  to  the  equation  of  state  must 
include  at  least  the  first  two  virial  coefficients.  Thermodynamic  data 
including  such  virial  corrections  are  usually  given  in  tabular  form 
(e.g.  references  (3)  -  (7))  or  in  the  graphical  form  of  a  Mollier 
diagram  (e.g.  reference  (81).  These  tabulations  and  diagrams  were  used 
to  provide  the  properties  of  nitrogen  at  high  densities  required  for 
calculations  concerning  tunnel  components  upstream  of  the  nozzle. 

Although  such  tabular  and  graphical  forms  are  very  useful  for  calcula¬ 
tions  done  by  hand,  they  are  not  very  well  suited  to  computer  applications 
such  as  the  computation  of  nozzle  contours. 

The  need  for  nitrogen  properties  at  moderately  high  densities  in 
a  form  suitable  for  computer  applications  was  filled  by  adapting  a 
versatile  computer  program  developed  by  Dr.  Harold  W.  Woolley  (reference 
(9))*.  This  program  can  compute  the  equilibrium  thermodynamic  properties 
of  arbitrary  mixtures  of  various  components**  at  moderately-high 
densities  by  treating  the  mixture  as  a  Cragoe-type  gas  with 

■jt - 

After  this  phase  of  the  work  was  completed,  Enkenhus  and  Culotta  (ref¬ 
erence  (10} )  presented  analytical  expressions  for  the  thermodynamic 
properties  of  nitrogen  which  agree  with  the  tabulations  of  Brahinsky 
(reference  (3))  to  better  than  one  percent. 

**Subroutines  for  the  six  components,  nitrogen,  oxygen,  argon,  neon, 
carbon  dioxide,  water  vapor,  were  included. 


2 


NOLTR  71-6 


In  Z 
P 


B'  +  C' p 


(1) 


where  Z  is  the  compressibility  of  the  gas  and  p  the  density.  The 
coefficients  B'  and  C'  have  been  determined  as  functions  of  tempera¬ 
ture  by  fitting  experimental  data.  Because  of  its  form,  this  expression 
leads  to  nitrogen  properties  which  agree  quite  well  with  published 
tabulations  (see  Figure  1)  provided  that  the  density  is  not  too  high. 

For  the  enthalpy  and  entropy  required  in  the  nozzle  supply  reservoir 
(i.e.  the  heater  vessel),  values  of  temperature,  pressure,  density, 
compressibility,  and  speed  of  sound  obtained  using  Woolley's  computer 
program  deviated  from  values  obtained  from  interpolation  of  Brahinsky's 
tabulations  (reference  (3))  by  less  than  0.4  percent,  1.2  percent,  and 
1.6  percent  for  the  Mach  number  10,  15,  and  20  nozzles,  respectively. 
Accordingly,  the  computer  program  was  judged  sufficiently  accurate 
and  was  used  to  provide  the  nitrogen  properties  required  for  computing 
the  nozzle  contours.  It  must  be  kept  in  mind  that  the  nozzle  calcula¬ 
tions  are  for  moderately-high  density  conditions.  For  higher  densities, 
such  as  in  the  gas  supply  vessels  upstream  of  the  heater,  Woolley's 
program  can  be  significantly  in  error  and  should  not  be  used  except 
with  the  utmost  caution. 

Given  values  for  any  two  thermodynamic  variables,  Woolley's 
program  will  compute  the  corresponding  values  for  most  other  thermo¬ 
dynamic  variables  of  interest.  Since  temperature  and  density  are  the 
natural  variables  for  the  equations  used  in  Woolley's  program,  a 
built-in  iterative  procedure  beginning  with  initial  guesses  for  tem¬ 
perature  and  density  is  used  for  other  pairs  of  variables.  The 
program  is  set  up  to  perform  multiple  calculations  either  by  reading 
successively  pairs  of  values  for  the  two  chosen  variables  or  by 
working  through  a  specified  range  of  the  two  variables  with  increments 
given  by  the  user. 

CONDENSATION  LINE  PROPERTIES 

In  order  to  attain  the  highest  test-section  Reynolds  number  for  a 
given  supply  pressure  and  test-section  Mach  number,  the  test-section 
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temperature  should  be  as  low  as  possible.  To  see  this,  we  anticipate 
the  results  somewhat  and  assume  that  the  temperature  and  pressure  in 
the  test  section  are  such  that  the  perfect  gas  law  applies,  the 
specific  heats  are  constant,  and  the  viscosity  is  a  function  of  tem¬ 
perature  only.  Then,  for  a  given  pressure  and  test-section  Mach 
number,  the  Reynolds  number  is  proportional  to  about  the  -1.25  power 
of  the  test-section  temperature  and,  therefore,  increases  as  the 
temperature  decreases.  This  behavior  is  evident  in  Figure  2  for  three 
supply  pressures  and  a  Mach  number  of  20. 

In  practice,  a  lower  limit  exists  for  the  test-section  temperature 
for  any  given  supply  pressure  because  the  gaseous  nitrogen  will 
condense  into  the  liquid  phase  if  cooled  sufficiently.  This  lower 
limit  is  the  temperature  at  which  the  nitrogen  is  on  the  verge  of 
condensing,  i.e.  the  condensation  temperature.  Operation  at  this 
limit  then  becomes  "condensation  threshold  operation."'  Using  the 
equilibrium  saturation  temperature-pressure  relation  for  the  onset  of 
condensation  given  by  Din  (reference  (4) ,  p.  86) 

log10  p  (mmHg)  =  6.49594  -  ^^^-6.6  {2) 

the  Reynolds  number  per  foot  for  condensation-threshold  operation  at 
Mach  number  20  is  shown  as  a  function  of  test-section  temperature  in 
Figure  2,  together  with  the  curves  for  constant  pressure  operation. 

The  figure  shows,  as  expected,  that  the  minimum  supply  pressure 
necessary  to  produce  a  given  Reynolds  number  per  foot  in  the  test 
section  is  the  pressure  corresponding  to  condensation  threshold  opera¬ 
tion. 

The  expression  given  by  Din  for  the  onset  of  condensation  applies 
only  if  the  nitrogen  does  not  become  supersaturated.  Daum  (reference 
(11))  shows  that  a  substantial  amount  of  supersaturation  can  occur 
when  a  test  gas  is  expanded  rapidly  to  high  velocity  and  low  pressure. 
If  the  distance  from  the  upstream  end  of  the  uniform  flow  region  to 
the  test  section  is  sufficiently  short,  the  test  gas  will  still  be  in 
a  supersaturated  state  when  it  reaches  the  test  section.  Applying 
the  procedure  outlined  by  Daum,  it  was  estimated  that  there  could  be 
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no  supexsaturation  in  the  Mach,  number  10  and  15  nozzles  and  only  a 
minimal  amount  in  the  Mach  number  20  nozzle.  Later  calculations 
based  on  the  final  nozzle  designs  agreed  with  the  early  estimates. 

For  the  Mach  number  20  case,  the  analysis  indicated  that  the 
test-section  static  temperature  could  be  lowered  about  2°K  before  the 
onset  of  condensation.  This  2°K  reduction  leads  to  a  five  percent 
decrease  in  the  required  supply  pressure  and  temperature  for  the  same 
Reynolds  number.  However,  it  was  decided  that  this  five  percent 
decrease  in  supply  conditions  was  not  sufficient  to  justify  risking 
the  possibility  of  condensation  in  the  test  section.  Moreover,  as 
will  be  seen  in  a  later  section,  the  Mach  20  nozzle  will  probably 
require  film  or  transpiration  cooling  to  protect  the  nozzle  throat 
during  the  test  run.  Such  cooling  may  reduce  the  stream  temperature 
and  precipitate  condensation.  This  possibility  is  reduced  somewhat  by 
operating  at  the  slightly  higher  temperature. 

Since  little  or  no  supersaturation  is  indicated,  Din's  expression 
was  chosen  to  represent  the  temperature-pressure  relation  at  condensa¬ 
tion  threshold.  If  a  temperature  is  given,  the  corresponding  pressure 
for  condensation  threshold  operation  is  obtained  from  Din's  expression. 
The  remaining  thermodynamic  properties  can  then  be  found  from  a 
Mollier  diagram,  perfect  gas  relations,  or  a  computer  program  such  as 
the  one  obtained  from  Woolley.  Woolley's  program  was  used  so  that  the 
condensation  line  properties  would  be  consistent  with  the  properties 
obtained  by  using  Woolley's  program  to  compute  isentropic  expansions 
from  the  selected  nozzle  supply  conditions.  The  pressure,  density, 
entropy,  enthalpy,  and  sound  speed  for  condensation  threshold  operation 
are  shown  in  Figures  3-7  and  listed  in  Table  1.  The  mass  flow  per  unit 
area  and  Reynolds  number  per  foot  are  shown  in  Figures  8  and  9 . 


VISCOSITY 


Low  Temperature.  The  viscosity  of  nitrogen  at  low  temperature 
required  for  Reynolds  number  calculations  was  obtained  from  the  dilute- 
gas  relation  given  by  Hirschf elder ,  Curtiss,  and  Bird  (reference  (12)) 


266.93X10~7  /MT  ,  g  . 
a2  fl  C2'2)*  CT*)  cm  sec 


C3) 
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where  T*  =  kT/e. 

For  nitrogen  M  =  28.016 

a  =  3.68 1A 
e/k  =  91. 46°K 


Substituting  these  values  and  changing  units 


y  =  7.0066X10 


-  7  /T 


nt2'2i"(T*) 


,  lbm  , 
'f tsec 


(4) 


■  •  (2  2 )  * 

The  collision  integral  0,  '  was  obtained  from  the  expression 


(T*)  =|  [6.771-6.38(T*-.5)+6.6(T*-.5)2-4(T*-.5)3] 


(5) 


which  represents  the  data  tabulated  in  reference  (12)  for  the  Lennard- 
Jones  (6-12)  potential  to  better  than  0.1  percent  accuracy  in  the 
range  30-80°K.  The  viscosity  calculated  from  these  expressions  (shown 
in  Figure  10  and  listed  in  Table  1)  was  used  in  the  Reynolds  number 
calculations  mentioned  previously. 

Higher  Temperature.  The  viscosity  of  nitrogen  at  higher  tempera¬ 
tures  and  densities  was  obtained  from  the  results  of  Brebach  and  Thodos 
(reference  (13))  who  give  a  reduced  state  correlation  for  the  viscosity 
of  diatomic  gases  at  atmospheric  pressure  and  a  density-dependent 
correction  to  account  for  higher  pressures.  Experimental  data  for  the 
viscosity  at  atmospheric  pressure  (y*)  are  correlated  by  the  expressions 


T. 


0.979 


R 


(Tr  <  1.0) 


(6) 


1.196 


0.659 


■R 


(Tr  k  3.5) 


(7) 


where  TR  =  T/Tc  and  Tc  is  the  critical  temperature.  For  nitrogen,  Tc 

is  given  as  126.2  degrees  Kelvin  and  y*T  as  .00865  centipoise.  Since 

c 
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Brebach  and  Thodos  do  not  give  an  expression  for  y*  in  the  range  1.0  < 
TR  <3.5,  a  Sutherland-type  formula  was  fitted  to  the  above  expres¬ 
sions  at  the  points  TR  =  1.0  and  TR  =  3.5  with  the  result 


rn  x  • 

1.7884  XR 
0.7884  +  Tr 


( 1 . 0<Tr<3 . 5) 


(8) 


The  viscosity  obtained  from  this  relation  seems  to  fit  the  data  shown 
by  Brebach  and  Thodos  to  better  than  0.1  percent  in  the  range  1.0  1 
Tr  <  3.5.  The  viscosity  obtained  using  the  appropriate  expressions 
given  above  for  the  three  temperature  ranges  is  shown  in  Figure  11. 

High  Density.  Pressures  higher  than  atmospheric  are  accounted 
for  by  a  density-dependent  correction  to  the  viscosity  at  atmospheric 
pressure  (y*)  as  given  by  the  above  expressions.  This  correction, 
obtained  by  fitting  the  nitrogen  data  shown  by  Brebach  and  Thodos  to 
an  accuracy  of  a  few  percent,  is  given  by 

y  -  y*  =  .01198  p  +  .06410p2  ( . 01^p£ . 35g/cm3)  (9) 

3 

where  y  and  y*  are  expressed  in  centipoises  and  p  in  g/cm  .  The  upper 

limit  on  the  density  range  is  dictated  by  the  accuracy  of  the  fit  but 

is  sufficient  for  present  purposes*.  The  lower  limit  coincides  with 

that  of  the  data  given  but  is  sufficient  also.  For  example,  in  the 

3 

nozzle  design  calculations,  the  density  of  .01g/cm  occurs  at  tempera¬ 
tures  greater  than  370°K.  At  this  temperature  and  density,  the 
viscosity  correction  is  0.6  percent  and  can  be  neglected. 

The  density-dependent  correction  to  the  viscosity  discussed  in 
the  previous  paragraph  is  non-negative.  Thus,  it  does  not  account 
for  the  decrease  in  viscosity  for  pressures  much  less  than  atmospheric 
as  shown  in  Figures  7  or  9  of  reference  (13).  However,  in  the  present 


*The  upper  limit  on  equation  (9)  exceeds  the  highest  supply  density 
(Mach  20  case)  by  30  percent.  For  cold  wall  nozzles,  the  local  density 
at  the  wall  can  exceed  the  supply  density.  For  the  present  nozzles, 
this  probably  does  not  occur,  but  at  worst  may  occur  only  in  the  very 
beginning  of  the  inlet  section  where  its  effects  can  be  neglected. 
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calculations,  the  temperature-pressure  combinations  occurring  are 
such  that  no  correction  of  this  type  is  needed. 

Application  of  Formulas.  The  viscosities  required  for  the  tur¬ 
bulent  boundary  layer  calculations  were  obtained  using  equations  (6) , 

(7) ,  and  (8)  for  the  viscosity  at  atmospheric  pressure  (y*) .  For 

3  3 

densities  less  than  .01g/cm  (.0194  slug/ft  ),  the  required  viscosity 

•  3 

is  equal  to  y*.  For  densities  greater  than  . Olg/cm  ,  the  density- 

dependent  correction  given  by  equation  (9)  is  also  used.  For  the  sake 

of  continuity,  the  viscosity  at  the  lower  temperatures  is  obtained  from 

the  Brebach-Thodos  expression  (equation  (6) )  rather  than  from  the 

Lennard-Jones  (6-12)  expression  (equation  (3)  or  (4) )  which  yields 

viscosities  about  2  percent  lower  than  the  Brebach-Thodos  expression 

over  the  range  30-80°K.  Thus,  the  Reynolds  numbers  near  the  nozzle 

exit  obtained  from  the  boundary  layer  calculation  will  be  about  2 

percent  lower  than  those  expected  on  the  basis  of  Figure  9. 


NOZZLE  OPERATING  CONDITIONS 


SUPPLY  CONDITIONS  FOR  CONDENSATION  THRESHOLD  OPERATION 

Given  the  conditions  along  Din's  condensation  line,  the  correspond¬ 
ing  nozzle  supply  (heater  vessel)  conditions  can  be  calculated  using 
Woolley's  computer  program.  Since  the  expansion  process  in  the  nozzle 
is  assumed  to  be  isentropic,  the  supply  entropy  is  given  by  Figure  5 
or  Table  1.  Then,  for  a  particular  condensation-line  temperature,  the 
supply  enthalpy  hQ  can  be  calculated  from 


h  2  2 

o  _  h  M  a 

R  R  2R 


(10) 


using  the  desired  Mach  number  and  the  condensation-line  enthalpy  and 
sound  speed  obtained  either  from  Figures  6  and  7  or  from  Table  1.  For 
each  supply  enthalpy-entropy  pair,  the  remaining  supply  conditions  can 
now  be  calculated  using  Woolley's  computer  program.  The  resulting 
supply  pressure  variations  with  condensation-line  temperature  are  shown 
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in  Figure  12  for  selected  Mach.  numbers  between  IQ  and  20.  The 
corresponding  supply  temperatures  and  densities  are  shown  in  Figures 
13  and  14,  respectively.  For  convenience,  in  Figure  15  the  supply 
pressure  is  shown  in  terms  of  the  Reynolds  number  per  foot  at  the 
nozzle  exit. 

CHOICE  OF  NOZZLE  SUPPLY  AND  TEST-SECTION  CONDITIONS 

The  particular  choice  of  supply  and  test  conditions  represents  a 
compromise  between  a  number  of  desirable  but  conflicting  features. 
Consideration  was  given  to  such  factors  as  test  duration,  Reynolds 
number  capability,  heater  size  and  strength,  nozzle  size,  pumping  and 
vacuum  requirements,  and  others,  in  order  to  find  the  optimum  operating 
conditions  which  could  be  attained  within  the  limitations  of  the 
allocated  funds.  Figures  16,  17,  and  18,  relating  the  supply  pressure, 
Reynolds  number,  mass  flow  rate,  volume  flow  rate  in  the  supply 
vessel,  and  the  nozzle  core  diameter,  were  made  to  simplify  the  search 
for  the  optimum  operating  conditions  at  Mach  number  10,  15,  and  20, 
respectively.*  The  nozzle  supply  and  test  section  conditions  resulting 
from  this  compromise  are  listed  in  Table  2. 

ACCURACY  OF  NITROGEN  PROPERTIES  AT  SUPPLY  CONDITIONS 

The  thermodynamic  properties  of  nitrogen  obtained  from  Woolley's 
computer  program  become  increasingly  inaccurate  at  very  high  densities 
because  only  two  coefficients  are  used  in  the  compressibility  equation 
(equation  (1) ) .  Figure  1  showed  that,  given  the  supply  pressure  and 
temperature,  the  calculated  density  is  quite  accurate.  However,  for 
the  nozzle  design  computations,  the  properties  must  be  calculated  at 
constant  entropy.  To  determine  the  accuracy  of  the  properties  for 
such  computations,  a  more  complete  comparison  was  made  with  the  tabulated 
results  of  Brahinsky  and  Neel  (reference  (3)).  To  avoid  the  need  for 
interpolations  in  both  the  entropy  and  temperature,  the  comparison  was 
made  at  the  design  entropy  and  at  the  temperature  in  Brahinsky' s  tab¬ 
ulations  nearest  to  the  design  supply  temperature.  The  results  of  this 

*Similar  figures  for  these  and  other  Mach  numbers  can  be  prepared  from 
the  data  given  in  Figures  8,  9,  12,  and  14. 
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comparison,  given  in  the  table  below  ,  show  that  the  nitrogen  properties 
obtained  from  Woolley's  program  are  sufficiently  accurate  for  the 
nozzle  design  computations.  Positive  percent  differences  indicate  that 
Woolley's  values  are  greater  than  those  obtained  from  Brahinsky  and 
Neel. 


Mach  Number 
Entropy  (S/R) 

Temperature  (supply)  °K 
Temperature  (comparison)  °K 
Pressure  difference  % 

Density  difference  % 

Enthalpy  difference  % 

Sound  Speed  difference  % 
Compressibility  difference  % 


10 

15 

20 

21.44 

22.15 

23.74 

1030 

1843 

2000 

1000 

1800 

2800 

-0.39 

-0.11 

-0.59 

1 

o 

■ 

to 

00 

-1.16 

-1.52 

-0.16 

-0.28 

0.05 

1 

o 

• 

o 

<1 

0.59 

1.10 

-0.17 

0.10 

0.92 

INVISCID  CORE  CALCULATION 

ISENTROPIC  EXPANSION  FROM  SUPPLY  CONDITIONS 

At  the  chosen  nozzle  supply  conditions,  the  properties  of  the 
nitrogen  test  gas  will  differ  substantially  from  the  properties  of  an 
ideal  gas  because  of  pressure  and  temperature  effects.  Thus,  the 
"real  gas"  properties  of  nitrogen  must  be  taken  into  account  in  the 
calculation  of  the  nozzle  contours.  Since  the  nitrogen  is  assumed  to 
remain  in  local  thermodynamic  equilibrium,  the  expansion  process  is 
isentropic.  The  nitrogen  properties'  can  then  be  calculated  and  tabulated 
before  beginning  the  nozzle  contour  calculations,  greatly  simplifying 
both  calculations. 

The  nitrogen  properties  for  the  isentropic  expansion  from  the 
chosen  supply  conditions  to  the  test  section  conditions  were  calculated 
using  Woolley's  equilibrium  thermodynamic  properties  computer  program 
mentioned  earlier.  Choosing  entropy  and  temperature  as  the  independent 
variables  for  Woolley's  program,  the  necessary  nitrogen  properties 
were  calculated  using  the  supply  entropy  together  with  temperatures 
ranging  from  the  supply  temperature  to  just  below  the  test  section 


10 


NOLTR  71-6 


temperature.  Since  a  constant  temperature  difference  becomes  a 
larger  percentage  difference  at  the  lower  temperatures,  the  temperature 
interval  was  decreased  as  the  test  section  temperature  was  approached. 
The  temperature  intervals  used  in  various  temperature  ranges  are 
tabulated  in  Table  3  for  the  three  nozzles. 

Woolley's  computer  program  was  modified  slightly  so  that  the 
desired  nitrogen  properties  were  punched  into  cards.  These  cards  were 
then  read  into  a  small  computer  program  which  transforms  Woolley's 
properties  and  writes  the  results  onto  magnetic  tape  in  the  form 
required  by  the  method  of  characteristics  program  used  for  the  inviscid 
core  calculations  in  the  supersonic  flow  region. 

SELECTION  OF  FLOW  REGIONS  FOR  CALCULATION  PURPOSES 

The  inviscid  flow  in  the  throat  region  of  a  supersonic  nozzle  can 
be  divided  into  the  four  regions  indicated  in  the  following  sketch. 
Region  I  represents  the  subsonic  flow  region  upstream  of  the  curved 
sonic  line  AA'.  Region  IV  represents  the  supersonic  flow  region 
downstream  of  a  characteristic  line  CC ' 


which  starts  an  arbitrary,  small  distance  downstream  of  the  sonic  point 
on  the  nozzle  centerline  and  proceeds  further  downstream  as  the  dis¬ 
tance  from  the  centerline  increases.  The  supersonic  flow  region 
upstream  of  the  characteristic  CC'  is  divided  into  Regions  II  and  III 
by  some  characteristic  BB *  of  the  family  opposite  to  that  of  CC ' . 

The  characteristic  BB '  represents  a  limit  upstream  of  which  the  method 
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of  characteristics  cannot  be  used  effectively.  Thus,  the  flow  in 
Region  III  can  be  calculated  by  the  method  of  characteristics. 

The  flow  in  Region  II  must  be  calculated  by  some  other  method, 
such  as  a  series  expansion  method.  Since  such  methods  have  not  yet 
been  developed  for  high- temperature ,  high-pressure  gases,  the  flow  in 
Region  II  is  not  calculated  but  treated  in  an  approximate  way  as 
discussed  in  the  subsonic  flow  section  following.  Although  the  flow 
in  Region  III  can  be  calculated  by  the  method  of  characteristics,  the 
region  is  so  small  that  it  is  treated  in  the  same  approximate  way  as 
Region  II.  For  example,  for  the  Mach  number  10  and  15  nozzles,  the 
point  C'  is  less  than  one-quarter  inch  downstream  of  the  sonic  point 
A.  For  the  Mach  number  20  nozzle,  the  distance  is  less  than  one-tenth 
of  an  inch.  For  the  present  nozzles,  there  is  little  likelihood  that 
the  contours  could  be  machined  accurately  enough  to  make  the  more 
exact  calculation  worthwhile. 

For  practical  purposes,  the  inviscid  flow  through  the  nozzle  can 
now  be  treated  in  two  parts:  The  supersonic  flow  region  downstream  of 
some  characteristic  CC '  and  the  combined  subsonic  flow  and  throat 
regions  upstream  of  CC 1 . 

SUPERSONIC  FLOW  REGION 

A  method  of  characteristics  procedure  for  calculating  the  inviscid 
core  contour  in  the  supersonic  flow  region  downstream  of  some  charac¬ 
teristic  such  as  CC '  Csee  sketch,  page  11)  for  uniform  flow  nozzles 
was  presented  in  reference  (14) .  The  equations  given  are  applicable 
to  the  isentropic  flow  of  real  and  ideal  gases  through  two-dimensional 
or  axisymmetric  nozzles.  A  FORTRAN  IV  computer  program  applyina  this 
procedure  to  the  flow  of  air  through  axisymmetric  nozzles  was  presented 
in  reference  (15)  under  the  name  Axisymmetric  Core  Program  (ACP) .  The 
ACP  program  is  written  in  such  a  way  that  the  gas  properties  enter  the 
program  in  the  form  of  tabulated  data  read  from  an  input  tape.  Thus, 
the  program  can  be  used  unchanged  for  any  gas  for  which  an  isentropic 
expansion  has  been  calculated  and  stored  on  tape  in  the  proper  form. 

The  utility  program  mentioned  in  a  previous  section  was  used  to  trans¬ 
form  the  isentropic  flow  data  obtained  from  Woolley's  program  and 
store  the  results  on  tape  in  the  form  required  for  the  ACP  program. 
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The  procedure  incorporated  into  the  ACP  program  requires  the 
specification  of  the  Mach  number  distribution  along  the  nozzle  center- 
line  from  the  sonic  point  to  the  beginning  of  the  uniform  flow  region 
near  the  nozzle  exit.  Initial  calculations  showed  that  the  two  center- 
line  Mach  number  distribution  formulas  included  in  the  program  as  given 
in  reference  (15)  were  unsuitable  for  use  in  the  Mach  number  15  and  20 
nozzles.  For  forty-foot  long  nozzles  at  these  Mach  numbers,  the  included 
Mach  number  formulas  lead  to  a  large  radius  of  curvature  at  the  nozzle 
throat  if  the  maximum  expansion  angle  is  kept  under  12  degrees.  Since 
the  heat  transfer  rates  near  the  throat  are  very  high,  such  long 
narrow  throat  regions  present  a  considerable  design  problem. 

To  decrease  the  throat  radius  of  curvature  and  consequently  the 
throat  heating  problem,  a  third  centerline  Mach  number  distribution 
formula  was  added  to  the  ACP  program.  This  formula  is  actually  a 
composite  using  the  Mach  number  distribution  corresponding  to  radial 
flow  with  two  patching  distributions,  one  from  the  sonic  point  to  the 
start  of  the  radial  flow  region,  and  one  from  the  radial  flow  region 
to  the  beginning  of  the  uniform  flow  region*.  Use  of  this  Mach  number 
distribution  allowed  the  radius  of  curvature  to  be  decreased  signifi¬ 
cantly  to  5.4  inches  (6  throat  diameters)  and  2.6  inches  (8  throat 
diameters)  for  the  Mach  number  15  and  20  nozzles,  respectively. 

Although  this  increased  the  peak  heat  transfer  at  the  nozzle  throat, 
it  reduced  the  total  heat  transfer  to  the  throat  region  and  decreased 
the  length  of  the  region  that  had  to  be  cooled.  Moreover,  the  nozzle 
throat  region  now  became  easier  to  machine  to  the  desired  accuracy. 

SUBSONIC  FLOW  AND  THROAT  REGION 

The  calculation  of  the  flow  in  the  subsonic  portion  of  the  nozzle 
for  a  high-temperature,  high-pressure  gas  is  very  difficult  at  best. 

The  same  is  true  for  much  of  the  small  region  of  supersonic  flow  in 
the  throat  region  (Regions  II  and  III  of  sketch,  page  11)  not  included 
in  the  method  of  characteristics  calculation  mentioned  in  the  last 

*Details  of  this  third  centerline  Mach  number  distribution,  its 
derivation,  instructions  for  its  use,  and  the  changes  made  to  ACP 
program  are  presented  in  Appendix  A. 
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section.  Consequently,  the  flow  in  the  subsonic  flow  end  throat 
regions  has  been  approximated  by  using  a  one-dimensional  flow  analysis, 
A  fourth-degree  polynomial  in  x  (axial  distance  from  sonic  point)  is 
used  to  define  the  inviscid  core  contour  in  these  regions.  Previously 
a  polynomial  was  used  to  define  the  nozzle  wall  contour  in  these 
regions.  The  present  approach  eliminates  the  sensitivity  of  the 
boundary  layer  calculation  in  the  supersonic  flow  region  to  the  choice 
of  initial  parameters.  Moreover,  although  the  present  approach  does 
not  necessarily  lead  to  a  more  accurate  wall  contour,  the  inviscid 
core  contour,  the  boundary  layer,  and,  thus,  the  nozzle  wall  contour 
are  each  on  a  firmer  theoretical  basis  than  previously. 

The  five  constants  in  the  fourth-degree  polynomial  are  determined 
by  specifying  the  length  of  the  subsonic  inlet  and  the  radius  and 
slope  at  the  upstream  end  of  the  inlet,  then  matching  at  the  throat 
(x  =  0)  the  radius,  slope  (  =  0) ,  and  curvature  obtained  from  the 
characteristics  calculation.  The  inlet  is  shown  schematically  in 
the  following  sketch. 


The  inviscid  core  radius  at  the  throat  r*  is  obtained  from  one-dimen¬ 
sional  mass  flow  considerations  and  is  supplied  by  the  ACP  characteristics 
calculation.  The  curvature  (actually  the  second  derivative  r*")  at 
the  throat  is  obtained  from  a  second-degree  polynomial  (quadratic) 
fitted  with  zero  slope  at  the  throat  to  the  throat  radius  r*  and  the 
radius  r^  of  the  characteristics  contour  at  some  point  x^.  The  inlet 
length  xq  (xq<0)  and  its  upstream  radius  rQ  and  slope  rQ  1  (  =  tan  oiQ) 
are  chosen  to  suit  the  particular  application.  For  the  present 
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nozzles,  the  radius  r  and  slope  r  1  were  chosen,  then  the  inlet 

o  o 

length  xq  was  varied  until  the  maximum  slope  of  the  inlet  was  about 
-30°.  The  inviscid  core  contour  in  the  subsonic  inlet  is  found  from 
the  polynomial 

r  =  Ax4  +  Bx3  +  Cx2  +  Dx  +  E  (11) 


where 
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The  inlet  area  ratio  is  given  by 


A 

A* 


(16) 


Cl  7) 


Since  the  flow  in  the  subsonic  inlet  and  the  throat  region  is  assumed 
to  be  one-dimensional,  flow  variables  at  various  locations  can  be 
obtained  from  the  isentropic  expansion  calculated  earlier,  once  the 
area  ratio  is  known.  However,  it  is  then  necessary  to  interpolate  six 
different  flow  variables  (pressure,  density,  enthalpy,  compressibility, 
sound  speed  and  Mach  number)  at  each  location  in  the  inlet.  Unless 
this  is  done  very  carefully,  the  consistency  between  the  flow  variables 
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will  be  lost  at  the  higher  area  ratios  where  the  variables  approach 
their  supply  values  asymptotically.  To  avoid  this  problem,  it  was 
decided  to  invert  the  polynomial  and  find  the  location  corresponding 
to  each  area  ratio  tabulated  in  the  isentropic  expansion.  The  inversion 
of  the  polynomial  was  carried  out  using  the  Newton-Raphson  iteration 
technique  and  proceeding  from  the  throat  upstream  to  the  beginning  of 
the  inlet. 

The  inviscid  core  contour  in  the  supersonic  portion  of  the  nozzle 
is  obtained  in  dimensionless  form  by  the  method  of  characteristics 
program  and  can  be  scaled  to  any  exit  diameter  desired.  This  scaling 
is  usually  done  in  the  boundary  layer  program.  The  calculation  of  the 
core  contour  for  the  subsonic  flow  and  throat  region  is  best  done  in 
dimensional  form  and  consequently  has  been  incorporated  into  the 
boundary  layer  program.  This  eliminates  the  need  for  a  separate 
program  and  simplifies  the  data  input  into  the  boundary  layer  program. 
Since  the  subsonic  portion  of  the  isentropic  expansion  table  input  into 
the  characteristics  program  does  not  end  up  on  the  output  tape,  the 
subsonic  data  is  read  into  the  boundary  layer  program  on  cards.  A 
program  which  punches  the  subsonic  portion  of  the  isentropic  expansion 
tape  onto  cards  in  the  proper  form  is  given  in  Appendix  B.  Selected 
coordinates  of  the  inviscid  core  contour  for  the  three  nozzles  are 
included  in  Table  4. 

TURBULENT  BOUNDARY  LAYER  CALCULATION 
GENERAL  THEORETICAL  APPROACH 

After  the  contour  of  the  inviscid  core  has  been  determined,  it 
must  be  corrected  in  order  to  account  for  the  growth  of  the  boundary 
layer  on  the  nozzle  wall.  This  correction  was  made  by  adding  to  the 
inviscid  core  contour  the  displacement  thickness  determined  using  the 
real  gas  turbulent  boundary  layer  calculation  method  presented  by 
Tetervin  in  reference  (16) .  In  addition  to  an  empirical  skin  friction 
law  and  an  integral  momentum  equation\alid  for  thick  as  well  as  thin 
boundary  layers,  Tetervin 's  method  utilizes  an  integral  moment-of- 
momentum  equation  to  calculate  the  momentum  thickness,  skin  friction, 
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and  a  parameter  which  determines  the  shape  of  the  velocity  profile. 

From  these  quantities  the  required  displacement  and  boundary  layer 
thicknesses  are  obtained.  Vector  addition  of  the  displacement  thick¬ 
ness  to  the  inviscid  core  contour  then  gives  the  contour  of  the  nozzle 
wall . 

A  FORTRAN  IV  computer  program  was  developed  to  apply  Tetervin's 
turbulent  boundary  layer  method  to  the  present  nozzle  design  calcula¬ 
tions.  The  program  as  developed  can  be  used  for  air  or  nitrogen  and 
is  compatible  with  the  method  of  characteristics  program  of  reference 
(15)  which  was  used  to  compute  the  supersonic  portion  of  the  inviscid 
core  contour.  For  convenience,  the  calculation  of  the  subsonic  portion 
of  the  inviscid  core  contour  described  in  the  preceding  section  was 
incorporated  into  the  boundary  layer  program.  The  user  has  the  option 
to  omit  the  inlet  core  calculation  and  begin  the  boundary  layer 
calculations  at  or  downstream  of  the  nozzle  throat.  The  program  and 
directions  for  its  use  are  given  in  Appendix  C.  In  the  present  cal¬ 
culations,  the  subsonic  inviscid  core  option  was  utilized.  Then  the 
turbulent  boundary  layer  growth  was  computed  from  the  beginning  of  the 
subsonic  inlet  to  the  exit  of  the  nozzle  using  the  empirical  Spalding- 
Chi  skin  friction  law  (reference  (17) )  and  the  Brebach-Thodos  viscosity 
correlation  for  nitrogen  discussed  previously. 

BASIC  EQUATIONS  AND  MODIFICATIONS 

The  basic  equations  used  are  Tetervin's  equations  (38)  to  (41), 

(45) ,  (51)  to  (54) ,  (59) ,  and  (61) .  With  the  exception  of  the 

equations  for  Aq,  Bq,  and  dH^/dx,  the  equations  are  programmed  as 
Tetervin  gives  them.  Using  Tetervin's  notation  and  defining 


(18) 


(19) 
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A  and  13  can  be  written  as 


(20) 


(21) 


Tetervin  concluded  that  the  most  reasonable  variation  of  H.  with 

l 


x  resulted  when  the  shear  stress  ratio 


t/t  was  taken  equal  to  unity 

W 


across  the  entire  boundary  layer.  Moreover,  the  most  plausible 

2 

results  were  obtained  when  the  first  (H .  -1)  in  the  equation  for 

2  2  1 

H.  )  with  H.  assumed  equal  to  its 
1f  1f 

Combining  these  results  and  defining 


dH./dx  was  replaced  by  (H,‘ 

i 
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(22) 


the  equation  for  dH^/dx  for  the  nonporous  wall  (v  =  0)  case  is 
obtained  in  the  form 


dH. 
_ l 

dx 


F  (x)  EL  (Hi+1) 


2 


+  H.  £, 
l  1 


-(ILtD  E2 


(23) 


In  the  derivation  of  his  equations,  Tetervin  uses  the  coordinates 

x  and  y  to  designate  distances  along  and  normal  to  the  wall  upon  which 

the  boundary  layer  is  growing.  For  the  nozzle  calculations,  it  is 

more  convenient  and  desirable  to  use  distances  measured  along  the 

nozzle  centerline  rather  than  along  the  nozzle  wall.  The  derivative 

with  respect  to  distance  s  (Tetervin1 s  x)  along  the  nozzle  wall  is 

w 

related  to  the  derivative  with  respect  to  distance  x  along  the  nozzle 
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centerline  by  the  cosine  of  the  wall  angle  a,  i.e. 


d  _  dsw  d _  1  d 

dx  ~  dx  ds  ~  cos  a  ds 

w  w 


(24) 


Thus,  Tetervin's  differential  equations  for  the  momentum  thickness  ©r 
and  the  transformed  shape  factor  (in  the  modified  form  given  in 
the  preceding  paragraph)  are  obtained  in  tegms  of  the  desired  center- 
line  derivatives  by  multiplying  through  by  . 


AUXILIARY  RELATIONS 

Wall  Temperature.  The  turbulent  boundary  layer  on  the  nozzle 
wall  can  be  significantly  affected  by  heat  transfer  to  or  from  the 
nozzle  wall.  To  uncouple  the  boundary  layer  development  from  the 
in-depth  temperature  response  of  the  nozzle  wall,  eliminating  the 
need  for  detailed  knowledge  of  the  wall  construction  and  material 
properties,  the  wall  surface  temperature  distribution  is  prescribed. 
Upstream  of  the  nozzle  throat,  the  wall  temperature  is  taken  to  be 
constant,  equal  to  the  throat  wall  temperature.  Downstream  of  the 
throat  (x  =  0) ,  the  temperature  distribution  is  given  by 


T  ( x ) 
w 


Ae-(Ex) 


c 


+  D 


(25) 


The  constants  in  this  relation  are  determined  by  specifying  the  tem¬ 
perature  of  the  nozzle  wall  at  the  nozzle  throat,  at  the  nozzle  exit, 
and  at  two  intermediate  points.  For  the  present  nozzles,  the  throat 
wall  temperature  was  taken  to  be  approximately  three- fourths  of  the 
supply  temperature.  The  exit  wall  temperature  was  taken  to  be  300°K, 
approximately  room  temperature.  The  intermediate  points  and  tempera¬ 
tures  chosen  are  indicated  in  Figure  19  which  shows  the  resulting  wall 
temperature  distributions  for  the  three  nozzles.  Although  this  pro¬ 
cedure  is  rather  arbitrary,  it  should  be  kept  in  mind  that,  especially 
in  the  throat  region  and  upstream,  the  nozzle  wall  which  was  initially 
near  room  temperature  will  be  heated  to  nearly  adiabatic  wall  con¬ 
ditions  (zero  heat  transfer).  Thus,  for  short  test  durations,  the 
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processes  are  essentially  time-dependent,  so  any  steady-state  wall 

temperature  calculation  is  also  rather  arbitrary.  The  form  of 

equation  (25)  for  T  is  based  on  the  assumption  that,  because  of 

w 

convective  heating,  the  wall  temperature  is  greatest  at  the  nozzle 
throat  and  decreases  continuously  toward  the  nozzle  exit.  Other  options 
available  in  the  computer  program  as  given  in  Appendix  C  are  for  a 
constant  wall  temperature  throughout  the  nozzle  and  for  adiabatic  wall 
conditions . 

Skin  Friction.  No  single  skin  friction  law  seems  to  give  good 
results  for  the  wide  range  of  flow  conditions  that  can  occur  in 
hypersonic  nozzles.  The  Ludwieg-Tillmann  law  (reference  (18)) 
includes  the  effect  of  profile  shape,  but  for  compressible  flows  must 
be  used  with  the  reference  enthalpy  method  (reference  (19)).  However, 
as  Tetervin  points  out  in  his  discussion,  the  reference  enthalpy 
method  can  give  skin  friction  coefficients  that  are  much  too  large  for 
very  cool  walls,  i.e.,for  walls  for  which  the  gas  enthalpy  at  the  wall 
(h  )  is  much  less  than  the  adiabatic  wall  enthalpy  (h  ) .  On  the 
other  hand,  the  5palding-Chi  skin  friction  law  does  not  account  for  the 
effect  of  profile  shape,  but  does  seem  to  give  fairly  accurate  results 
when  h ^/h  is  small.  For  the  present  calculations,  the  dependence  on 

hw/haw  was  considered  the  more  important,  so  the  Spalding-Chi  law  has 
been  used.  However,  as  these  two  skin  friction  laws  tend  to  complement 
each  other,  both  have  been  incorporated  into  the  computer  program. 

Prandtl  Number.  The  variation  of  the  Prandtl  number  is  difficult 
to  calculate  for  most  gases  and  conditions  of  interest  in  nozzle  design. 
Consequently,  a  constant  value  of  Prandtl  number  Pr  equal  to  .72  has 
been  used  in  the  present  calculations  and  incorporated  into  the  com¬ 
puter  program.  The  use  of  a  constant  Prandtl  number  results  in 
constant  values  for  such  Prandtl  number  related  quantities  as  the 
recovery  factor  and  the  Reynolds  analogy  factor.  Thus, 

h  -  h  1/3 

recovery  factor  =  a^_— ^ --  =  Pr  =  .  89628  (26) 

o  e 

where  h  is  the  gas  enthalpy  and  the  subscripts  aw,  e,  and  o  refer  to 
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adiabatic  wall,  edge  of  the  boundary  layer,  and  supply  conditions, 
respectively.  Also,  following  Colburn, 

Reynolds  analogy  factor  =  '/2  ~  Pr  =  1.244  8  C27) 

where  St  is  the  Stanton  number  and  the  skin  friction  coefficient. 

GAS  PROPERTIES 

The  gas  properties  necessary  for  the  boundary  layer  calculation 
are  made  available  to  the  computer  program  in  three  ways.  The  supply 
and  throat  conditions  and  most  of  the  necessary  conditions  along  the 
edge  of  the  boundary  layer  for  the  portion  of  the  nozzle  downstream 
of  the  throat  are  input  from  the  tape  produced  by  the  Axisymmetric 
Core  Program  (ACP)  discussed  previously.  The  corresponding  conditions 
for  the  nozzle  portion  upstream  of  the  throat  are  input  from  cards  as 
described  in  Appendix  C.  In  addition  to  these  constant  entropy 
properties,  the  density  distribution  through  the  boundary  layer,  the 
viscosity  at  the  edge  of  the  boundary  layer,  and  the  enthalpy  (or  tem¬ 
perature  for  the  adiabatic  wall  case)  and  viscosity  at  the  wall  are 
needed  to  complete  the  calculations.  Thus,  expressions  for  calculating 
these  quantities  must  be  incorporated  into  the  computer  program. 

Since  the  available  expressions  often  are  not  in  terms  of  the 
desired  variables,  some  procedure,  usually  iterative,  must  be  used 
with  the  expressions  to  obtain  the  desired  results.  For  example,  the 
pressure  through  the  boundary  layer  is  assumed  to  be  constant,  equal 
to  the  known  value  pg  at  the  edge  of  the  boundary  layer.  The  enthalpy 
distribution  through  the  boundary  layer  is  calculated  from  the  Crocco 
relation  using  the  velocity  distribution.  Thus,  the  most  natural 
variables  for  calculating  the  required  density  distribution  are  the 
pressure  p  (=Pe)  and  the  enthalpy  h.  In  functional  form  this  may  be 
written  as  p (p,  h) .  If  an  expression  of  the  form  h(p,  p)  is  used  (as 
in  the  case  of  air) ,  then  an  iterative  procedure  must  be  used  which, 
from  an  initial  guess  for  the  density  p,  calculates  h,  and  corrects  p 
based  on  the  difference  between  the  desired  and  calculated  values  of 
h,  repeating  this  procedure  until  the  desired  degree  of  convergence  is 


21 


NOLTR  71-6 


attained.  If,  on  the  other  hand,  an  expression  of  the  form  hCT,  p)  is 
available  Cas  in  the  case  of  nitrogen) ,  then  a  different  iterative 
procedure,  this  one  involving  the  equation  of  state  p  (T,  p) ,  is  needed. 
Since  the  available  expressions  for  air  and  nitrogen  differ  and  require 
different  supporting  procedures  to  achieve  the  desired  results,  the 
two  cases  will  be  discussed  separately. 

Nitrogen.  It  was  pointed  out  in  the  preceding  paragraphs  that 
expressions  must  be  incorporated  into  the  computer  program  for  cal¬ 
culating  needed  gas  properties  not  input  from  cards  or  tape.  For 
nitrogen,  the  expressions  used  are 

p  =  pZRT  (28) 


h 

R 


(z  +  §>  T  + . ■ . 

exp  (jjj)  -1 


(29) 


Z 


1  (X  =  log, n  (-£-)  <  0) 

u  PST 

1  +  .00482X2eX  (0<X<2.3) 


(30) 


where  0  =  3390°K  and  p  =  2.423516x10  2  slugs/ft2.  The  expression 

O  -L 

for  the  compressibility  Z  was  developed  by  Harris  (reference  (20) ) 
using  the  thermodynamic  data  of  Lewis  and  Neel  (references  (21)  and 
(22))  and  is  valid  for  pressures  up  to  1000  atmospheres  and  temperatures 
up  to  3000°K.  For  most  applications  this  range  is  adequate.  The 
pressure  in  the  inlet  and  throat  regions  6f  the  Mach  20  nozzle  and  in 
part  of  the  Mach  15  inlet  exceeds  1000  atmospheres,  but  in  these 
regions  the  boundary  layer  is  thin  and  the  contour  somewhat  arbitrary. 
For  these  reasons,  no  attempt  was  made  to  extend  the  range  to  higher 
pressures . 

Given  the  pressure  p  (=p  )  and  the  enthalpy  h,  these  equations 
can  be  solved  using  the  Newton-Raphson  technique  to  obtain  the  density, 
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temperature,  and  compressibility.  This  procedure  is  used  to  calculate 
the  density  distribution  through  the  boundary  layer,  the  viscosities 
at  the  edge  of  the  boundary  layer  and  at  the  wall  (using  equations  (6) 
to  (9)),  and  the  wall  temperature  for  the  case  of  an  adiabatic  wall. 

The  same  equations  can  be  solved,  given  the  pressure  and  wall  tem¬ 
perature,  by  another,  slightly  different,  Newton-Raphson  procedure  to 
find  the  enthalpy  at  the  wall. 

Air .  Air  properties  required  for  the  boundary  layer  calculation, 
but  not  input  into  the  program  from  cards  or  tape,  must  be  calculated 
from  appropriate  expressions  incorporated  into  the  program.  To 
calculate  the  density  distribution  in  the  boundary  layer  from  the 
pressure  p  (=Pe)  and  enthalpy  h,  i.e.,  to  calculate  p(p,h),  a  form  of 
Newton's  rule  of  false  position  is  used  with  an  expression  for  enthalpy 
of  the  functional  form  h(p,p).  Beginning  with  an  initial  guess  for  p, 
the  procedure  is  repeated  until  the  desired  degree  of  convergence  is 
attained.  The  expression  for  h(p,p)  plus  another  for  the  compressibility 
Z  in  the  form  Z(p,p)  were  given  originally  by  Grabau  (reference  (23)) 
and  were  corrected,  evaluated,  and  presented  as  the  subroutines  ENTHLP 
and  COMP  in  reference  (24) .  The  use  of  these  enthalpy  and  compres¬ 
sibility  expressions  together  with  the  equation  of  state  allows  the 
wall  temperature  for  the  case  of  an  adiabatic  wall  to  be  calculated 
from  the  pressure  and  enthalpy. 

The  wall  enthalpy  could  be  calculated  from  the  pressure  and  tem¬ 
perature  using  these  expressions,  but  the  procedure  is  rather  involved. 
Since  the  supply  pressures  for  air  nozzles  are  moderate  and  the  wall 
temperatures  are  not  too  high,  the  enthalpy  is  assumed  to  be  a  function 
of  temperature  only  and  calculated  from 


h  3.3175xl0-4  (T-499.17) 

RT  J.  l-expl-. 012287  (T-499.17) ] 


(31) 


This  expression  is  a  fit  to  the  ideal  gas  properties  (pressure  effects 
neglected)  of  air  and  is  good  for  temperatures  in  degrees  Kelvin  up  to 
approximately  2000°K.  The  remaining  quantities  needed,  the  viscosities 
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at  the  wall  and  at  the  edge  of  the  boundary  layer,  are  calculated  from 
fits  of  the  form  y(p,h)  to  the  viscosity  data  of  Hansen  (reference 
(25) )  . 


RESULTS  AND  DISCUSSION 


The  more  important  results  of  the  calculations  of  the  nitrogen 
properties,  the  isentropic  expansion,  and  the  inviscid  core  contour 
have  been  presented  in  earlier  sections  of  this  report.  In  reference 
(16) ,  Tetervin  presented  and  discussed  in  detail  most  of  the  results 
of  the  boundary  layer  calculations  for  the  three  nozzles.  Therefore, 
in, this  section,  only  the  remaining  aspects  of  the  boundary  layer 
calculation  and  the  final  nozzle  contours  are  considered. 

INITIAL  INPUT  CHOICES 

Before  the  boundary  layer  calculations  can  be  carried  out,  several 
choices  must  be  made  from  the  options  available.  For  the  computer 
program  as  given  in  Appendix  Cf  the  user  must  decide  whether  to  use  the 
Spalding-Chi  or  the  Ludwieg-Tillmann  skin  friction  law,  whether  to  use 
adiabatic  wall  conditions  or  a  prescribed  wall  temperature  distribution 
whether  to  include  the  calculation  of  the  subsonic  inlet,  and  which 
numerical  method  to  use  to  integrate  the  differential  equations.  The 
user  must  then  specify  certain  input  data  appropriate  to  the  choice  of 
options  being  used,  e.g.,the  wall  temperature  at  four  locations  if  the 
variable  wall  temperature  distribution  is  being  used,  or  the  length  and 
initial  radius  and  slope  of  the  subsonic  inlet,  if  it  is  being  included 
For  the  present  nozzle  calculations,  the  Spalding-Chi  skin  friction 
law  was  used,  the  variable  wall  temperature  distributions  shown  in 
Figure  19  were  prescribed,  and  the  subsonic  inlet  calculations  were 
included. 

The  numerical  integration  of  the  differential  equations  for  the 
momentum  thickness  6^  and  the  shape  factor  EL  was  performed  using  the 
utility  subroutine  FN0L2  (reference  (26) )  in  the  predictor-corrector 
mode.  To  begin  the  integration,  values  of  0  and  II  ^  had  to  be 
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specified  at  the  initial  station.  These  values  can  be  specified 
arbitrarily.  Alternatively,  can  be  obtained  by  setting  it  equal  to 
the  flat  plate  value  calculated  from  Tetervin's  equation  (61) 


using  the  specified  value  of 


This  procedure  eliminates  the  need 


for  specifying  the  initial  value  of  II.  and  was  used  in  the  present 
calculations . 

If  the  boundary  layer  calculation  is  started  at  the  nozzle  throat 
position,  the  choice  of  the  initial  values  for  0^  and  affect  the 
calculated  boundary  layer  thicknesses  downstream  of  the  throat,  possibly 
all  the  way  to  the  nozzle  exit.  However,  by  beginning  the  calculation 
upstream  of  the  throat,  this  effect  can  be  reduced.  In  particular, 
this  effect  can  be  minimized  by  beginning  the  calculation  at  the 
start  of  the  subsonic  inlet.  This  result  is  the  major  justification 
for  including  the  approximate  treatment  of  the  subsonic  inlet  described 
earlier.  For  each  of  the  present  nozzles,  the  subsonic  inlet  cal¬ 
culation  was  included  to  minimize  the  effect  of  the  initial  value  of 

0  .  Figure  20  shows  the  variation  of  0  in  the  inlet  for  initial 
r  r 

values  differing  by  two  orders  of  magnitude  for  the  Mach  number  15 
nozzle.  The  curves  converge  to  a  single  minimum  value  just  upstream 
of  the  nozzle  throat.  The  same  behavior  is  observed  for  Ih.  Thus, 
the  boundary  layer  calculation  downstream  of  the  nozzle  throat  is 
independent  of  the  choice  of  initial  value  for  0^  within  a  rather  wide 
range.  Moreover,  in  the  upstream  end  of  the  inlet,  the  boundary  layer 
is  very  thin  compared  to  the  inlet  radius,  so  that  the  changes  in  the 
nozzle  wall  contour  in  this  region  are  negligibly  small.  Thus,  the 
calculated  wall  contour  for  the  entire  nozzle  is  unaffected  by  the 
initial  value  chosen  for  0  . 

Values  of  the  various  parameters  and  initial  values  used  in  the 
final  calculations  are  given  in  Table  5. 


PRINCIPAL  DIMENSIONS  OF  NOZZLES 

In  order  to  be  able  to  use  the  same  test  cell  and  model  support 
system  for  all  three  nozzles,  the  three  nozzles  were  designed  to  have 
the  same  exit  diameters  and  the  same  lengths.  Preliminary  calculations 
indicated  that  a  throat-to-exit  length  of  about  45  feet  and  an  exit 
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diameter  of  about  5  feet  would  be  necessary  for  the  Mach  number  20 
nozzle  if  the  maximum  expansion  half-angle  was  to  be  kept  near  12  degrees. 
It  was  decided  to  design  the  nozzles  such  that  at  a  distance  of  40  feet 
from  the  nozzle  throat  the  nozzle  diameter  would  be  five  feet.  The 
remaining  portion  of  the  nozzle  would  not  be  fabricated,  since,  if 
this  cut-off  portion  was  not  too  long,  the  nozzle  performance  would 
not  be  appreciably  affected.  The  optimum  cut-off  point  is  that  axial 
location  at  which  the  boundary  layer  enters  the  uniform  flow  region 
inside  the  Mach  cone  (see  Figure  21) .  Downstream  of  this  point,  the 
usable  test  core  decreases,  even  though  the  cross-sectional  area  of 
the  nozzle  continues  to  increase.  Iterations,  involving  the  boundary 
layer  program  only,  were  sufficient  to  produce  contours  for  the  Mach 
number  10  and  15  nozzles  in  which  the  optimum  cut-off  point  occurred 
within  a  few  inches  of  the  40- foot  station.  These  small  differences 
made  further  iterations  unnecessary  and  the  nozzles  were  cut  off  at 
the  40- foot  station.  The  length  of  the  cut-off  portion  was  1.72  and 
1.50  feet  for  the  Mach  number  10  and  15  nozzles,  respectively.  For 
the  Mach  number  20  nozzle,  the  optimum  cut-off  point  occurred  at  the 
44.7-foot  station  or  1.88  feet  from  the  end  of  the  uncut  nozzle.  To 
move  the  optimum  cut-off  point  nearer  to  the  40-foot  station  would 
require  a  maximum  expansion  angle  appreciably  greater  than  any  considered 
acceptable.  For  this  reason,  no  further  iterations  were  made  for  the 
Mach  number  20  nozzle.  The  diameter  of  the  usable  test  core  was  18.2 
inches  at  the  40- foot  cut-off  point  rather  than  the  21  inches  which 
would  have  been  attained  if  the  nozzle  had  been  cut  off  at  the  optimum 
location.  The  seven  square  foot  test  core  thus  obtained  was  considered 
adequate.  Principal  dimensions  of  the  three  nozzles  are  given  in 
Table  6. 

HEAT  TRANSFER  TO  NOZZLE  WALL 

The  boundary  layer  calculations  mentioned  above  assume  that  the 
nozzle  wall  temperature  can  be  kept  to  a  reasonable  level  without 
transpiration  or  film  cooling.  This  can  be  done  for  the  Mach  number 
10  nozzle  and  probably  for  the  Mach  number  15  nozzle.  However,  heat 
transfer  calculations  for  the  Mach  number  20  nozzle  show  that  massive 


26 


NQLTR  71-6 


transpiration  or  film  cooling  may  be  necessary  in  the  throat  region 

to  keep  the  nozzle  throat  from  eroding  early  in  the  test  run. 

The  nozzle  wall  temperature  distributions  for  the  three  nozzles 

were  assumed  to  vary  from  about  three-fourths  of  the  supply  temperature 

at  the  throat  to  300°K  at  the  nozzle  exit  as  shown  in  Figure  19.  The 

convective  heat  transfer  from  the  boundary  layer  to  the  nozzle  wall 

was  calculated  from  the  skin  friction  coefficient  using  Colburn's 

form  of  Reynolds  analogy  as  given  by  equation  (27) .  The  throat  region 

is  of  particular  interest  since  the  heat  transfer  rate  is  much  higher 

there  than  elsewhere  in  the  nozzle.  For  example,  for  the  assumed 

temperature  distributions,  the  peak  heat  transfer  rate  occurs  near  the 

2 

nozzle  throat  and  is  770,  19600,  and  46000  BTU/ft  sec  for  the  Mach 

number  10,  15,  and  20  nozzles,  respectively,  while  at  the  nozzle  exit, 

2 

the  rate  is  about  1  BTU/ft  sec  for  each  nozzle.  The  heat  transfer 
rates  in  the  throat  regions  of  the  three  nozzles  are  shown  in  Figure 
22. 

For  the  Mach  number  20  nozzle,  the  heat  transfer  rate  is  extremely 
high  in  the  throat  region.  Coupling  this  with  the  rather  high  value 
of  1600°K  assumed  for  the  wall  temperature  and  the  high  static  pressure 
of  1560  atmospheres  at  the  throat,  it  is  clear  that  the  design  of  the 
throat  region  for  the  Mach  20  nozzle  presents  a  formidable  problem. 

A  number  of  possible  solutions  to  this  problem  are  now  being  studied. 

The  heat  transfer  rate  in  the  throat  region  of  the  Mach  number 
15  nozzle  is  down  to  one-half  the  rate  for  the  Mach  20  nozzle.  More¬ 
over,  the  wall  temperature  and  pressure  levels  are  more  moderate.  Thus, 
there  is  a  reasonable  chance  that  the  throat  section  of  the  Mach  15 
nozzle  can  be  made  to  have  an  acceptable  lifetime  without  using  trans¬ 
piration  or  film  cooling. 

Conditions  of  temperature,  pressure,  and  heat  transfer  rate  in  the 
throat  region  of  the  Mach  10  nozzle  are  sufficiently  moderate  that  no 
transpiration  or  film  cooling  is  needed. 

NOZZLE  WALL  CONTOURS 

The  boundary  layer  program  discussed  in  this  report  does  not 
include  provisions  for  using  transpiration  or  film  cooling  to  protect 
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the  nozzle  during  operation.  In  the  previous  section,  it  was  shown 
that  the  Mach  number  10  and  15  nozzles  can  probably  operate  acceptably 
without  such  cooling,  but  most  likely  the  Mach  number  20  nozzle  can 
not.  Since  calculations  which  take  such  cooling  into  account  have  not 
yet  been  completed,  the  Mach  20  nozzle  has  been  calculated  with  the 
present  program  as  though  no  such  cooling  is  required.  Thus,  the 
Mach  20  nozzle  contour  is  considered  preliminary  at  the  present  time. 

During  the  boundary  layer  .calculation,  the  inviscid  core  contour 
(subsonic  and  supersonic  flow  regions)  is  automatically  corrected  to 
account  for  the  effect  of  the  boundary  layer  and  the  resulting  nozzle 
contour  is  tabulated.  However,  inspection  of  the  calculated  contours 
showed  that  each  of  the  three  nozzle  contours  required  some  smoothing 
in  a  small  portion  of  the  subsonic  inlet.  The  reason  for  this  has  not 
yet  been  determined.  The  corrected  nozzle  wall  radii,  obtained  by 
hand-smoothing  the  calculated  results,  are  given  in  Table  7  for  each  of 
the  nozzles.  Selected  coordinates  of  the  nozzle  contour,  together 
with  the  coordinates  of  the  inviscid  core,  are  presented  in  Table  4. 

After  the  Mach  number  10  and  15  nozzle  contours  had  been  submitted 
for  fabrication,  an  error  was  found  in  the  boundary  layer  computer 
program.  The  density-dependent  correction  to  the  viscosity  was  being 
calculated  incorrectly.  The  program  was  corrected  and  the  contours 
recalculated.  The  nozzle  wall  was  found  to  be  displaced  outward 
slightly  with  the  increase  in  radius  ranging  from  zero  at  the  throat 
to  about  one-thousandth  of  an  inch  at  the  nozzle  exit.  Since  these 
changes  were  quite  small  compared  to  the  uncertainties  inherent  in  the 
boundary  layer  calculation,  the  changes  were  neglected.  Thus,  the 
nozzle  wall  coordinates  given  in  Table  4  are  the  uncorrected  values 
used  in  making  the  nozzles.  The  computer  program  given  in  Appendix  C 
is  the  corrected  version. 

ACCURACY  OF  CONTOURS 

The  nozzle  design  method  presented  in  this  report  is  subject  to 
error  from  a  number  of  sources. 

1.  The  assumption  of  local  thermodynamic  equilibrium  and 

isentropic  flow  in  the  inviscid  core  is  only  an  approximation. 
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For  the  present  nozzle  design  calculations,  this  approximation 
should  be  very  good.  However,  for  more  severe  flow  conditions 
and  extreme  nozzle  shapes,  this  approximation  can  lead  to  errors 
which  are  no  longer  negligibly  small. 

2.  Inaccuracies  in  the  gas  properties  for  the  flow  conditions 
occurring  in  the  nozzle  can  affect  the  calculated  results  sub¬ 
stantially.  Again,  for  the  present  nozzles,  errors  from  this 
source  are  considered  small. 

3.  The  assumptions  of  one-dimensional  flow  in  the  subsonic  inlet 
and  throat  regions  can  cause  errors  which  are  very  difficult  to 
evaluate.  For  perfect  gas  flows,  better  (but  more  complicated) 
methods  can  be  used.  However,  these  methods  are  not  yet  applicable 
to  flows  in  which  high  temperature  or  high  pressure  effects  on 

the  gas  properties  are  important. 

4.  Errors  introduced  by  the  method  of  characteristics  procedure 
used  for  the  supersonic  inviscid  core  and  by  the  numerical 
procedures  used  in  the  boundary  layer  program  can  be  reduced  to 

a  negligible  level  by  using  sufficiently  small  step  sizes.  There 
is  no  difficulty  here. 

5.  The  main  sources  of  error  in  the  boundary  layer  calculation 
seem  to  be  in  the  accuracy  of  the  available  skin  friction  laws 
for  highly  cooled  walls  and  in  the  assumed  wall  temperature 
distribution.  For  the  design  of  nozzles  (such  as  the  present 
ones)  having  thick  boundary  layers,  the  errors  in  the  boundary 
layer  calculation  probably  outweigh  all  others.  Fortunately,  as 
Tetervin  points  out,  errors  in  the  boundary  layer  calculation 
are  at  least  partially  self -correcting.  If,  for  example,  the 
boundary  layer  is  actually  thicker  than  calculated,  the  flow  Mach 
number  will  be  less  than  the  design  value.  However,  for  the  same 
nozzle  geometry  and  supply  conditions,  a  decrease  in  the  Mach 
number  leads  to  a  decrease  in  the  boundary  layer  thickness  which 
in  turn  acts  to  increase  the  Mach  number.  Thus,  any  tendency  of 
the  boundary  layer  to  thicken  or  thin,  brings  about  a  change  in 
the  Mach  number  distribution  which  opposes  that  tendency.  This 
equilibrium-like  behavior  tends  to  reduce  the  effect  of  inaccu¬ 
racies  in  the  boundary  layer  calculation. 
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6.  The  final  source  of  error  in  the  nozzle  contour  is  in  the 
actual  fabrication  of  the  nozzle.  Prohibitive  costs  and  the 
capabilities  of  available  equipment  combine  to  limit  the  accuracy 
with  which  the  nozzle  can  be  made. 

In  summary,  the  final  nozzle  is  subject  to  errors  from  a  number 
of  sources.  Some  of  these  errors  can  be  made  negligibly  small,  but 
others  cannot,  at  least  at  the  present  time.  Because  of  the  number 
and  nature  of  the  possible  errors,  it  is  not  possible  to  obtain  a 
numerical  estimate  of  the  overall  accuracy  of  the  calculated  nozzle 
contours.  However,  the  three  nozzle  contours  presented  and  discussed 
in  this  report  are  considered  to  be  as  accurate  as  is  now  possible  to 
make  them . 


CONCLUDING  REMARKS 


The  procedure  used  in  designing  axisymmetric  nozzles  for  the  NOL 
Hypervelocity  Wind  Tunnel  and  the  resulting  nozzles  for  operation  at 
Mach  numbers  10,  15,  and  20  with  high-density  nitrogen  and  thick 
turbulent  boundary  layers  have  been  presented  in  detail. 

The  NOL  Hypervelocity  Wind  Tunnel  is  a  unique  aerodynamic  facility 
for  hypersonic  research  and  developmental  testing.  The  capability  of 
achieving  a  high  Mach  number  -  high  Reynolds  number  flow  with 
relatively  long  run  times  and  large  model  sizes  will  permit  the 
accumulation  of  large  amounts  of  aerodynamic  test  data,  greatly 
facilitating  the  design  and  development  of  advanced,  high-performance 
re-entry  vehicles. 

This  report  summarizing  the  nozzle  design  for  the  facility  is 
one  of  a  series.  Other  technical  reports  documenting  various  aspects 
of  the  aerodynamic  and  mechanical  design,  tunnel  performance,  and 
test  results  obtained  will  be  published  as  the  material  becomes 
available. 
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TABLE  1 

NITROGEN  PROPERTIES  AT  CONDENSATION  THRESHOLD 


SOUND 


TEMPERATURE 

PRESSURE 

DENSITY 

ENTROPY* 

ENTHALPY 

SPEED 

0  K 

ATM 

AMAGAT 

BTU/LBM 

FT/SEC 

40.0 

9.0327E-05 

6 • 1656E-04 

25.312 

17.755 

422.96 

40.1 

9.5210E-05 

6.4826E-04 

25.269 

17.800 

423.4b 

40.2 

1 .0033E-04 

6.8140E-04 

25.225 

17.845 

424.01 

40,3 

1 .0568E-04 

7. 1600E-04 

25.182 

17.889 

424.54 

40.4 

1.1 129E-04 

7.52 1 4E-0  4 

25.139 

17.934 

425.07 

40.5 

1.1716E-04 

7.8986E-04 

25.096 

17.979 

425.59 

40.6 

1. 2331 E- 04 

8.2923E-04 

25.053 

18.023 

426.12 

40.7 

1  •  2973E-04 

8.7030E-04 

25.011 

18.068 

426.64 

40.8 

1 .3645E-04 

9.131 5  E-0  4 

24.969 

18.113 

427.16 

40.9 

1 .4348E-04 

9.5783E-04 

24.928 

18.158 

427.69 

41.0 

1 . 5082E-04 

1 . 0044E-0  3 

24.886 

18.202 

428.21 

41.1 

1 . 58  50E-04 

1 .0529E-03 

24.845 

18.247 

428.73 

41.2 

1 . 6652F-04 

1.1035E-03 

24.804 

18.292 

429.25 

41.3 

1 .7489E-04 

1  ,  1  562  E-0  3 

24.764 

18.336 

429.77 

41.4 

1 .8363E-04 

1.211 1 E-0  3 

24.723 

18.381 

430.29 

41.5 

1 .9276E-04 

1.2682E-03 

24.683 

18.426 

430.81 

41.6 

2.0228E-04 

1.3  2  77E-0  3 

24.643 

18.470 

431.33 

41.7 

2.1222E-04 

1.3895E-03 

24.604 

18.515 

431.85 

41.8 

2. 2258 E- 04 

1.4539E-03 

24.565 

18.560 

432.36 

41.9 

2  •  3339E-04 

1 . 5209E-03 

24.526 

18.604 

432.88 

42.0 

2.4465E-04 

1.5905E-03 

24.487 

18.649 

433.39 

42.1 

2  «  5639E-04 

1 .6628E-03 

24.448 

18.694 

433.91 

42.2 

2.6P62E-04 

1 .7380E-03 

24.410 

18.738 

434.42 

42.3 

2. 8137 E- 04 

1 . 8 1 62  E-0  3 

24.372 

18.783 

434.94 

42.4 

2  •  9464E-04 

1 .8974E-03 

24.334 

18.828 

435.45 

42.5 

3  •  0846E— 04 

1 .9817E-03 

24.296 

18.872 

435.96 

42.6 

3  •  2284E-04 

2.0693E-03 

24.259 

18.917 

436.48 

42.7 

3.3781E-04 

2.1601E-03 

24.222 

18.961 

436.99 

42.8 

3, 5339 E- 04 

2.2545E-03 

24.185 

19.006 

437.50 

42.9 

3 . 69  59E-04 

2.3523E-03 

24.148 

19.051 

438.01 

43.0 

3 • 8644E-04 

2 • 4539E-0  3 

24.112 

19.095 

438.52 

43.1 

4.0396E-04 

2.5592E-03 

24.076 

19.140 

439.03 

43.2 

4.2217E-04 

2.6683E-03 

24.040 

19.185 

439.53 

43.3 

4.4109E-04 

2.7815E-03 

24.004 

19.229 

440.04 

43.4 

4.6076E-04 

2.8988E-03 

23.968 

19.274 

440.55 

43.5 

4.8118E-04 

3 • 0204E-0  3 

23.933 

19.319 

441.06 

43.6 

5.0240E-04 

3.1463E-03 

23.898 

19.363 

441.56 

43.7 

5 . 2443  E-04 

3 • 2  768  E-0  3 

23.863 

19.408 

442.07 

43.8 

5,4729 E-04 

3.4119E-03 

23.828 

19.453 

442.57 

43.9 

5. 7103 E-04 

3.5517E-03 

23.794 

19.497 

443.07 

44.0 

5. 9566 E— 04 

3.696  5  E-0  3 

23.760 

19.542 

443.5b 

44.1 

6.212 1 E-04 

3 • 8  464E-0  3 

23.726 

19.586 

444.08 

44.2 

6. 4771 E-04 

4.0014E-03 

23.692 

19.631 

444.58 

44.3 

6,75 19E-04 

4.16 1 8E-0  3 

23.658 

19.676 

445.08 

44.4 

7.0368E-04 

4.3277E-03 

23.625 

19.720 

445.59 

44.5 

7 . 3322E-04 

4 . 499  2  E-0  3 

23.591 

19.765 

446.09 

44.6 

7  •  6383E-04 

4.676  5  E-0  3 

23.558 

19.810 

446.59 

44.7 

7 ,9555 E-04 

4.8599E-03 

23.526 

19.854 

447.09 

44.8 

8. 2841 E-04 

5.0493E-03 

23.493 

19.899 

447.58 

44.9 

8  •  6244E-04 

5 • 2  45 1 E-0  3 

23.460 

19.943 

448.08 

VISCOSITY 
LBM/FT  SEC 

1.8470E-06 
1.8514E-06 
1 • 8  55  7E-0  6 
1.8601E-06 
1.8644E-06 
1.8688E-06 
1.8731E-06 
1.8775E-06 
1.8818E-06 
1.8862E-06 
1.8906E-06 
1. 8949E-06 
1.8993E-06 
1. 9037E-06 
1.9080E-06 
1 • 9 1 24E-06 
1.9168E-06 
1.9212E-06 
1. 9256E-06 
1.9299E-06 
1.9343E-06 
1.9387E-06 
1.9431E-06 
1.9475E-06 
1.9519E-06 
1.9563E-06 
1.9607E-06 
1.9651E-06 
1 • 9695E-06 
1.9739E-06 
1.9783E-06 
1.9827E-06 
1.9871E-06 
1.9915E-06 
1.9959E-06 
2.0004E-06 
2.0048E-06 
2.0092E-06 
2.0136E-06 
2.0180E-06 
2.0225E-06 
2. 0269E-06 
2.0313E-06 
2.0357E-06 
2.0402E-06 
2.0446E-C6 
2 • 0490  E-06 
2.0535E-06 
2.0579E-06 
2.0624E-06 


*  ENTROPY  S  EXPRESSED  IN  DIMENSIONLESS  FORM  S/R 
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TABLE  1  (CONT.) 


SOUND 


TEMPERATURE 

PRESSURE 

DENSITY 

ENTROPY* 

ENTHALPY 

SPEED 

°K 

ATM 

AMAGAT 

BTU/LBM 

FT/SEC 

45.0 

8.9769 E-04 

5.4473E-03 

23.428 

19.988 

448.58 

45.1 

9.3418E-04 

5.6562E-03 

23.396 

20.033 

449.08 

45.2 

9.7195E-04 

5.8719E-03 

23.364 

20.077 

449.57 

45.3 

1.0110E-03 

6. 0946 E-03 

23.332 

20.122 

450.07 

45.4 

1.0515 E-03 

6.3245E-03 

23.301 

20.166 

450.56 

45.5 

1 .0933E-03 

6.561 8  E-0  3 

23.270 

20.211 

451.06 

45.6 

1 . 1366E-03 

6.8067E-03 

23.238 

20.256 

451.55 

45.7 

1.1814E-03 

7.0  5  93  E-0  3 

23.207 

20.300 

452.04 

45.8 

1  •  22  77  E- 03 

7.31 99E-0  3 

23.177 

20.345 

452.54 

45.9 

1 .2755E-03 

7.5887E-03 

23.146 

20.389 

453.03 

46.0 

1.3250E-03 

7.8659E-03 

23.116 

20.434 

453.52 

46.1 

1 .3761E-03 

8. 1517E-03 

23.085 

20.479 

454.01 

46.2 

1 .4289E-03 

8 • 446  2  E-0  3 

23.055 

20.523 

454.50 

46.3 

1 .4835E-03 

8 . 749  8  E-0  3 

23.025 

20.568 

454.99 

46.4 

1 . 5398E-03 

9.0626E-03 

22.996 

20.612 

455.48 

46.5 

1.5980E-03 

9.3850E-03 

22.966 

20.657 

455.97 

46.6 

1.6581 E-03 

9.7170E-03 

22.937 

20.701 

456.46 

46.7 

1,7201 E-03 

1 .0059E-02 

22.907 

20.746 

456.94 

46.8 

1.7841 E-03 

1.0411E-02 

22.878 

20.791 

457.43 

46.9 

1 . 8502E-03 

1.0774E-02 

22.849 

20.835 

457.92 

47.0 

1 .9183E-03 

1.1 1 47E-0  2 

22.821 

20.880 

458.40 

47.1 

1 • 9887E-03 

1,15  3 1 E-0  2 

22.792 

20.924 

458.89 

47.2 

2.0612E-03 

1.1926E-02 

22.764 

20.969 

459.37 

47.3 

2.1  360E-03 

1 .2333E-02 

22.735 

21.013 

459.86 

47.4 

2.2131 E-03 

1 . 2752E-02 

22.707 

21.058 

460.34 

47.5 

2.2926E-03 

1 .3182E-02 

22.679 

21.102 

460.82 

47.6 

2  •  3746E-03 

1,362  5E-02 

22.652 

21.147 

461.30 

47.7 

2.4591E-03 

1.4080E-02 

22.624 

21.191 

461.79 

47.8 

2  •  5461E-03 

1 .4548E-02 

22.596 

21.236 

462.27 

47.9 

2 . 63  58  E-03 

1.5029E-02 

22.569 

21.280 

462.75 

48.0 

2. 7282 E-03 

1.5524E-02 

22.542 

21.325 

463.23 

48.1 

2.8234E-03 

1.6032E-02 

22.515 

21.369 

463.71 

48.2 

2.9214E-03 

1,655  5  E-0  2 

22,488 

21.414 

464.19 

48.3 

3 . 02 2  3E-03 

1 . 709 1 E-0  2 

22.461 

21.458 

464.66 

48.4 

3. 1262 E-03 

1.7642E-02 

22.435 

21.503 

465.14 

48.5 

3. 2331 E-03 

1.8208E-02 

22.408 

21.547 

465.62 

48.6 

3.3431E-03 

1 .8789E-02 

22.382 

21.592 

466.09 

48.7 

3 . 45  64E-03 

1 .9386E-02 

22.356 

21.636 

466.57 

48.8 

3 . 5729E-03 

1 .9O99E-02 

22.330 

21.681 

467.05 

48.9 

3. 6928 E-03 

2.0628E-02 

22.304 

21.725 

467.52 

49.0 

3. 8161 E-03 

2.1274E-02 

22.278 

21.770 

467.99 

49.1 

3  •  9429E-03 

2. 1936E-02 

22.253 

21.814 

468.47 

49.2 

4.0733E-03 

2.2616E-02 

22.227 

21.858 

468.94 

49.3 

4. 2073 E-03 

2.33 1 3E-0  2 

22.202 

21.903 

469.41 

49.4 

4.3451 E-03 

2 . 402  8  E-0  2 

22.177 

21.947 

469.88 

49.5 

4  «  4868  E-03 

2.476  2  E-02 

22.152 

21.992 

470.36 

49.6 

4.6324E-03 

2.5514E-02 

22.127 

22.036 

470.83 

49.7 

4  •  78  2  0E-03 

2 . 628  5  E-0  2 

22.102 

22.080 

471.30 

49.8 

4.9357E-03 

2.7076E-02 

22.078 

22.125 

471.77 

49.9 

5 .0936E-03 

2.78  8  7E-0  2 

22.053 

22.169 

472.23 

VISCOSITY 
LBM/FT  SEC 

2.0668 E-06 
2.0713E-06 
2.0757E-06 
2 • 0  802  E-06 
2.0846 E-06 
2.0891 E-06 
2.0935E-06 
2.0980E-06 
2.1 024E-06 
2. 1069E-06 
2. 1114E-06 
2. 1158E-06 
2. 1203E-06 
2. 1247E-06 
2. 1292E-06 
2. 1337E-06 
2.1 382E-06 
2.1 426E-06 
2.147 1 E-06 
2.1516E-06 
2. 1561 E-06 
2. 1605E-06 
2 . 1650E-06 
2 . 1695E-06 
2. 1740E-06 
2. 1785E-06 
2. 1830 E-06 
2. 1875E-06 
2.1 920E-06 
2. 1966 E-06 
2.2  0 10  E-06 
2.2054E-06 
2 . 2099E-06 
2.2 144E-06 
2.2190E-06 
2.2235E-06 
2.2280E-06 
2.2325E-06 
2.2370E-06 
2.241 5E-06 
2. 2460E-06 
2.2505E-06 
2. 2550E-06 
2 . 2595E-06 
2 . 2641 E-06 
2 . 2686E-06 
2.2731 E-06 
2 . 2776E-06 
2.2821E-06 
2 . 2867E-06 


*  ENTROPY  S  EXPRESSED  IN  DIMENSIONLESS  FORM  S/R 
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TABLE  1  (CONT.) 


SOUND 


TEMPERATURE 

PRESSURE 

DENSITY 

ENTROPY* 

ENTHALPY 

SPEED 

VISCOSITY 

0  K 

ATM 

AMAGAT 

BTU/LBM 

FT/SEC 

LBM/FT  SEC 

50.0 

5 

. 2558 E- 03 

2. 8718 E-02 

22.029 

22.214 

472.70 

2.291 2  E-06 

50.1 

5 

.4223E-03 

2.9569E-02 

22.004 

22.258 

473.17 

2.2957E-06 

50.2 

5 

• 5934E-03 

3 • 0442  E-0  2 

21.980 

22.302 

473.64 

2.3003E-06 

50.3 

5 

• 7690E-03 

3.1336E-02 

21.956 

22.347 

474.11 

2 . 3  048  E-0  6 

50.4 

5 

• 9493E-03 

3.2252E-02 

21.933 

22.391 

474.57 

2 . 3093E-06 

50.5 

6 

•1344E-03 

3. 3190 E-02 

21.909 

22.435 

475.04 

2.3138E-06 

50.6 

6 

• 3244E-03 

3.41 50E-0  2 

21.885 

22.480 

475.50 

2.3 1 84E-06 

50.7 

6 

•  5  I 93E-03 

3.5 1 34E-02 

21.862 

22.524 

475.97 

2 . 3229E-06 

50.8 

6 

•  71 93  E-03 

3.6 1 42  E,-0  2 

21.838 

22.568 

476.43 

2.3275E-06 

50.9 

6 

•  9245  E-03 

3.717  3  E-0  2 

21.815 

22.613 

476.89 

2.3320E-06 

51.0 

7 

•1350E-03 

3 • 8229E-02 

21.792 

22.657 

477.36 

2 . 3365E-06 

51.1 

7 

• 3509E-03 

3.9309E-02 

21.769 

22.701 

477.82 

2.341 1 E-06 

51.2 

7 

• 5724E-03 

4.0415E-02 

21.746 

22.746 

478.28 

2 • 3456E-06 

51.3 

7 

. 7995  E-03 

4.1 547E-02 

21.723 

22.790 

478.74 

2 . 3502E-06 

51.4 

8 

.0323 E-03 

4. 2705 E-02 

21.701 

22.834 

479.20 

2.3547E-06 

51.5 

8 

.2710E-03 

4.38  8  9E-0  2 

21.678 

22.878 

479.66 

2. 3593 E-06 

51.6 

8 

.51 57E-03 

4. 5101E-02 

21.656 

22.923 

480.12 

2.3638E-06 

51.7 

8 

.7665 E-03 

4. 6341 E-02 

21.634 

22.967 

480.58 

2.3684E-06 

51.8 

9 

.0235 E-03 

4,7608 E-02 

21.611 

23.011 

481.03 

2.3729E-06 

51.9 

9 

.2868 E-03 

4,8905  E-02 

21.589 

23.055 

481.49 

2 . 3775E-06 

52.0 

9 

.5567 E-03 

5.0230E-02 

21.567 

23.100 

481.95 

2.3820E-06 

52.1 

9 

.8331 E-03 

5.158  5  E-0  2 

21.546 

23.144 

482.41 

2.3866E-06 

52.2 

1 

,01 16E-02 

5 • 2970E-0  2 

21.524 

23.188 

482.86 

2.3911E-06 

52.3 

1 

. 0406E-02 

5.4386E-02 

21.502 

23.232 

483.32 

2.3957E-06 

52.4 

1 

. 0703E-02 

5.5833E-02 

21.481 

23.276 

483.77 

2.4002E-06 

52.5 

1 

, 1007 E- 02 

5,731 2E-02 

21.459 

23.320 

484.22 

2.4048E-06 

52.6 

1 

. 1319E-02 

5.8823E-02 

21.438 

23.365 

484.68 

2 • 4094E-06 

52.7 

1 

.1638E-02 

6. 0366 E-02 

21.417 

23.409 

485.13 

2.4139E-06 

52.8 

1 

• 1964E-02 

6. 1943 E-02 

21.396 

23.453 

485.58 

2.4185E-06 

52.9 

1 

.2298E-02 

6,35  5  4E-0  2 

21.375 

23.497 

486.03 

2.423 1 E-06 

53.0 

1 

• 2640E-02 

6.5199E-02 

21.354 

23.541 

486.48 

2 . 4276E-06 

53.1 

1 

• 2990E-02 

6.6880E-02 

21.333 

23.585 

486.94 

2 • 4322E-06 

53.2 

1 

• 3348E-02 

6.8596E-02 

21.313 

23.629 

487.38 

2 • 4368  E-06 

53.3 

1 

•3714E-02 

7 . 0348E-02 

21.292 

23.673 

487.83 

2.4413E-06 

53.4 

1 

• 4089E-02 

7.21 3  6E-0  2 

21.272 

23.717 

488.28 

2 . 4459  E-0  6 

53.5 

1 

. 4472  E-02 

7 .3962 E-02 

21.251 

23.762 

488.73 

2.4505E-06 

53.6 

1 

• 4864E-02 

7. 5826 E-02 

21.231 

23.806 

489.18 

2.4550E-06 

53.7 

1 

.5265 E-02 

7.772  8  E-0  2 

21.211 

23.850 

489.62 

2 . 4596E-06 

53,8 

1 

. 5674E-02 

7 , 9669E-02 

21.191 

23.894 

490.07 

2.4642E-06 

53.9 

1 

. 6094E-02 

8.165 1 E-0  2 

21.171 

23.938 

490.52 

2 • 4688E-06 

54.0 

1 

•6522E-02 

8.3672E-02 

21.151 

23.982 

490.96 

2.4733E-06 

54.1 

1 

•6960E-02 

8.5734E-02 

21.131 

24.026 

491.40 

2.4779E-06 

54.2 

1 

. 7408  E-02 

8.783  8  E-0  2 

21.112 

24.070 

491.85 

2.4825E-06 

54.3 

1 

. 7865E-02 

8 . 9984E-02 

21.692 

24.113 

492.29 

2. 4871 E-06 

54.4 

1 

• 8333E-02 

9,2172  E-0  2 

21.073 

24.157 

492.73 

2.4917E-06 

54.5 

1 

•8810E-02 

9 . 440  4E-0  2 

21.053 

24.201 

493.18 

2.4963E-06 

54.6 

1 

. 9299E-02 

9.6680E-02 

21.034 

24.245 

493.62 

2. 5008E-06 

54.7 

1 

.9797 E-02 

9.9001E-02 

21.015 

24.289 

494.06 

2 • 5054E-06 

54.8 

2 

.0307 E-02 

1.01 37  E-0 1 

20.996 

24.333 

494.50 

2.5 100L-06 

54.9 

2 

•0827E-02 

1 .0378E-01 

20.977 

24.377 

494.94 

2.5146E-06 

*  ENTROPY  S  EXPRESSED  IN  DIMENSIONLESS  FORM  S/R 
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TABLE  1  (CONT.) 


SOUND 


TEMPERATURE 

PRESSURE 

DENSITY 

ENTROPY* 

ENTHALPY 

SPEED 

°  K 

ATM 

AMAGAT 

BTU/LBM 

FT/SEC 

55.0 

2.1359E-02 

1 .0624E-01 

20.958 

24.421 

495.38 

55.1 

2 . 1901E-02 

1 .0874E-01 

20.939 

24.465 

495.81 

55.2 

2 . 2456E-02 

1 .1130E-01 

20.920 

24.508 

496.25 

55.3 

2.3021 E-02 

1 . 1390E-01 

20.902 

24.552 

496.69 

55.4 

2  •  3599E-02 

1.16  5  5  E-0 1 

20.883 

24.596 

497.13 

55.5 

2.4189E-02 

1 . 1926E-01 

20.865 

24.640 

497.56 

55.6 

2.4791E-02 

1.2201E-01 

20.846 

24.684 

498.00 

55.7 

2 . 5405F-02 

1.2482E-01 

20.828 

24.727 

498.43 

55.8 

2  •  6033E-02 

1 .2767E-01 

20.810 

24.771 

498.86 

55.9 

2  •  6672E-02 

1 .3058E-01 

20.792 

24.815 

499.30 

56.0 

2.732  5E-02 

1.335  5E-0 1 

20.774 

24.858 

499.73 

56.1 

2  •  7992E-02 

1 .3657E-01 

20.756 

24.902 

500.16 

56.2 

2  •  8671 E-02 

1 .3964E-01 

20.738 

24.946 

500.59 

56.3 

2 . 9365E-02 

1 .4277E-01 

20.720 

24.989 

501.03 

56.4 

3 . 00 72E-02 

1 .4596E-01 

20.702 

25.033 

501.46 

56.5 

3  •  0793E-02 

1 .4920E-01 

20.685 

25.077 

501.89 

56.6 

3.1  529E-02 

1.5  2  50E-0 1 

20.667 

25.120 

502.31 

56.7 

3 . 2279E-02 

1 .5586E-01 

20.650 

25.164 

502.74 

56.8 

3 . 3044E-02 

1 .5929E-01 

20.633 

25.207 

503.17 

56.9 

3  •  3P24E-02 

1.6277E-01 

20.615 

25.251 

503.60 

57,0 

3 . 46 1 9  E-02 

1 , 663 1 E-0 1 

20.598 

25.294 

504.02 

57.1 

3.5430E-02 

1.6991 E-0 1 

20.581 

25.338 

504.45 

57.2 

3.6256E-02 

1 . 7358E-0  1 

20.564 

25.381 

504.87 

57.3 

3. 7098 E-02 

1 .773 1 E-0 1 

20.547 

25.425 

505.30 

57.4 

3 . 7956E-02 

1,811 1 E-0 1 

20.530 

25.468 

505.72 

57.5 

3, 8831 E-02 

1 .8497E-01 

20.513 

25.512 

506.15 

57.6 

3 • 9722E-0? 

1 .8890E-01 

20.497 

25.555 

506.57 

57.7 

4.0630E-02 

1 .9289E-01 

20.480 

25.599 

506.99 

57.8 

4. 1555 E-02 

1 .9666E-01 

20.464 

25.642 

507.41 

57.9 

4, 2498E-02 

2.0109E-01 

20.447 

25.685 

507.83 

58.0 

4. 3458 E-02 

2.0529E-01 

20.431 

25.729 

508.25 

58 . 1 

4  •  4436  E-02 

2 • 0956E-0 I 

20.414 

25.772 

508.67 

58.2 

4. 5432 E-02 

2.1390E-01 

20.398 

25.815 

509.09 

58.3 

4 . 6446E-02 

2.183  2  E-0 1 

20.382 

25.859 

509.51 

58.4 

4. 7479 E-02 

2.2280E-01 

20.366 

25.902 

509.92 

58.5 

4. 8531 E-02 

2.273  7E-0 1 

20.350 

25.945 

510.34 

58.6 

4.9602E-02 

2.3200E-01 

20.334 

25.988 

510.76 

58.7 

5 . 069 2E-02 

2,367 1 E-0  1 

20.318 

26.031 

511.17 

58.8 

5.1802E-02 

2,41 50E-0  1 

20.302 

26.075 

511.59 

58.9 

5 . 2932E-02 

2.4637E-01 

20.286 

26.118 

512.00 

59.0 

5 . 408  2  E-02 

2.513 1 E-0 1 

20.271 

26.161 

512.42 

59.1 

5 • 5252E-02 

2.563  3  E-0 1 

20.255 

26.204 

512.83 

59.2 

5 • 6443E-02 

2.61 44E-0 1 

20.239 

26.247 

513.24 

59.3 

5  •  7656E-02 

2 • 666  2  E-0  1 

20.224 

26.290 

513.65 

59.4 

5  •  8889E-02 

2.718  8  E-0 1 

20.209 

26.333 

514.06 

59.5 

6.01 44E- 02 

2.7723E-01 

20.193 

26.376 

514.47 

59.6 

6.142  IE-02 

2 . 82  66E-0  1 

20.178 

26.419 

514.88 

59.7 

6.2720E-02 

2.881 8  E-0 1 

20.163 

26.462 

515.29 

59.8 

6.4042E-02 

2 . 9  379  E-0  1 

20.148 

26.505 

515.70 

59.9 

6. 5386 E-02 

2.9947E-01 

20.133 

26.548 

516.10 

VISCOSITY 
LBM/FT  SEC 

2.5192E-06 
2.5238E-06 
2.5284E-06 
2. 5330E-06 
2.5376E-06 
2.5421E-06 
2.5467E-06 
2.5513E-06 
2.5  559E-06 
2. 5605E-06 
2.5651 E-06 
2. 5697E-06 
2. 5743 E-06 
2. 5789E-06 
2.5835E-06 
2. 5881E-06 
2.5927E-06 
2 • 5974E-06 
2 • 6020E-06 
2.6066E-06 
2.6112E-06 
2. 6158 E-06 
2 • 6204E-06 
2.6250E-06 
2 . 6296E-06 
2.6342E-06 
2 . 6388E-06 
2 . 643  5  E-06 
2.6481E-06 
2.6527E-06 
2.6573E-06 
2 . 66 1 9E-06 
2 • 6666E-06 
2. 6712 E-06 
2.6758E-06 
2 . 6804E-06 
2 . 6850E-06 
2 . 6897E-06 
2 . 6943E-06 
2.6989E-06 
2.7035E-06 
2 . 7082E-06 
2. 7128 E-06 
2.71 74E-0  6 
2.7220E-06 
2 • 726  7E-06 
2.73 13E-06 
2.7359E-06 
2.7  406E-06 
2 • 745  2  E-06 


*  ENTROPY  S  EXPRESSED  IN  DIMENSIONLESS  FORM  S/R 
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NOZZLE  SUPPLY 

TABLE 

AND  TEST 

2 

SECTION 

CONDITIONS 

MACH  NUMBER 

10 

15 

20 

SUPPLY 

PRESSURE 

430. 

2365. 

3110. 

ATM 

TEMPERATURE 

1030. 

1843. 

2800. 

0  K 

1854. 

3317. 

5040. 

°R 

ENTROPY  (S/R) 

21.44 

22.15 

23.74 

ENTHALPY 

492.4 

1016. 

1594. 

BTU/LBM 

DENSITY 

7.721 

18.68 

16.92 

LBM/FT3 

TEST  SECTION 

PRESSURE  X  102 

1.127 

0  •  448 

.0612 

ATM 

TEMPERATURE 

52.58 

49.50 

44 .06 

0  K 

94.64 

89.10 

79.31 

°R 

DENSITY  X  103 

4.567 

1.929 

0.296 

LBM/FT3 

VELOCITY 

4846. 

7055. 

8877. 

ET/SEC 

REYNOLDS  NUMBER 

9.2 

6.0 

1.3 

PER  FOOT  X  10 
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TABLE  3 

INTERVALS  USED  WITH  WOOLLEY'S  PROGRAM 


EROM 

MACH  10  450  ATM 

1029.80  ° K 
427  ATM 
1020  ° K 
890  0  K 
840  0  K 
60.0  °K 
53.0  °K 
52.5  ° K 


MACH  15  2380  ATM 

1842.38  °K 
1840  0 K 
1600  °K 
1520  °K 
6  5.0  0  K 
50.0  °K 


MACH  20  3130  ATM 

2799.51  °K 
3106  ATM 
2797.6  °K 
2790  ° K 
2400  °K 
2320  °  K 
70.0  °K 
44.96  °K 


(INTERVAL)  TO 


(-1  ATM) 

430  ATM 

(-.05  °K) 

1028.50  °  K 

(-1  ATM) 

420  ATM 

(-5  °K) 

890  °K 

(-1  °K ) 

840  °K 

(-5  °K) 

60.  °K 

(-.25  °K) 

53.0  °K 

(-.02  °K) 

52.5  °K 

(-.25  0 K ) 

40.25  °K 

(-1  ATM) 

2365 

ATM 

(-.02  °K) 

1840. 

92  0  K 

(-5  °K) 

1600 

°  K 

(-1  °  K ) 

1520 

°  K 

V 

o 

iTi 

I 

65.0 

°  K 

(-.25  °K) 

50.0 

°  K 

(-.02  0 K ) 

49.2 

°  K 

-1  ATM) 

3110  ATM 

o 

CO 

o 

o 

1 

2798.76  °  K 

-1  ATM) 

3101  ATM 

-.3  *»K) 

2790.4  ° K 

(-5  °K) 

2400  ° K 

(-1  °K) 

2320  °K 

(-5  °K) 

70.0  °K 

-.25  °K) 

45.25  °K 

-.02  °K) 

43.5  ° K 
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DTSTANCF 
FROM 
THROAT 
<  FT) 

-.75 

-.70 

-.65 

-.60 

-.55 

-.50 

-.45 

-.40 

-.38 

-.36 

-.34 

-.32 

-.30 

-.28 

-.26 

-.24 

-.22 

-.20 

-.19 

-.18 

-.17 

-.16 

-.15 

-.14 

-.13 

-.12 

-.11 

-.10 

-.09 

-.08 

-.07 

-.06 

-.05 

-.04 

-.03 

-.02 

-.01 

0.00 

0.01 

0.02 

0.03 

0.04 


TABLE  4 

SELECTED  INVISCID  CORE  AND  NOZZLE  WALL  COORDINATES 


MACH  NUMBER  =  10 

MACH  NUMBER  =  15 

MACH  NUMBER  =  20 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

OF  WALL 

OF  CORE 

OF  WALL 

OF  CORE 

OF  WALL 

OF  CORE 

(FT) 

(  FT  ) 

(FT  ) 

(  FT  ) 

(FT) 

(FT) 

0.33958 

0.33924 

0.33112 

0.33077 

0.31348 

0.31317 

0.29083 

0.29057 

0.26512 

0.26490 

0.23781 

0.23764 

0.21092 

0.21078 

0.19167 

0.19140 

0.18543 

0.18531 

0.18130 

0.18106 

0.16667 

0.16657 

0.17567 

0.17406 

0.17385 

0.16416 

0.16403 

0.16633 

0.16572 

0.16554 

0.15970 

0.15955 

0.15772 

0.15761 

0.15651 

0.15635 

0.15329 

0.15314 

0.15002 

0.14991 

0.14666 

0.14652 

0.14533 

0.14519 

0.14285 

0.14275 

0.13639 

0.13627 

0.13639 

0.13626 

0.13619 

0.13608 

0.12592 

0.12582 

0.12660 

0.12648 

0.13011 

0.13000 

0.11545 

0.11536 

0.11600 

0.12463 

0.12452 

0.10492 

0.10517 

0.11975 

0.11964 

0.09467 

0.09417 

0.09409 

0.11548 

0.11537 

0.08508 

0.08306 

0.08299 

0.11357 

0.11346 

0.08067 

0.07758 

0.07751 

0.11181 

0.11170 

0.07650 

0.07217 

0.07211 

0.11019 

0.11008 

0.07225 

0.07219 

0.06685 

0.06680 

0.10872 

0.10860 

0.06842 

0.06837 

0.06166 

0.06161 

0.10738 

0 .10726 

0.06472 

0.06467 

0.05662 

0.05658 

0.10618 

0. 10606 

0.06124 

0.06119 

0.05171 

0.05167 

0.10511 

0.10498 

0.05796 

0.05791 

0.04706 

0.04702 

0.10416 

0.10403 

0.05489 

0.05484 

0.04261 

0.04257 

0.10333 

0.10320 

0.05205 

0.05200 

0.03838 

0.03835 

0.10262 

0.10249 

0.04944 

0.04940 

0.03443 

0.03440 

0.10202 

0.10188 

0.04708 

0.04704 

0.03075 

0.03072 

0.10151 

0.10137 

0.04497 

0.04493 

0.02680 

0.02678 

0. 101 1  0 

0.1O095 

0.04311 

0.04307 

0.02429 

0.02426 

0.10078 

0.10062 

0.04151 

0.04146 

0.02161 

0.02159 

0.10052 

0.10036 

0.04016 

0.04011 

0.01928 

0.01926 

0.10036 

0.10020 

0.03906 

0.03901 

0.01733 

0.01731 

0.10024 

0 .10007 

0.03821 

0.03816 

0.01578 

0.01576 

0.10017 

0.09999 

0.03761 

0.03756 

0.01465 

0.01463 

0.10014 

0.09995 

0.03726 

0.03721 

0.01396 

0.01394 

0.10014 

0 . 09994 

0.03715 

0.03709 

0.01373 

0.01370 

0.10016 

0.00995 

0.03727 

0.03721 

0.01395 

0.01392 

0.10019 

0.09998 

0.03761 

0 .03754 

0.01441 

0.01436 

0. "10023 

0.10001 

0.03804 

0.03796 

0.01517 

0.01511 

0.10029 

0.10006 

0.03859 

0.03850 

0.01607 

0.01600 
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TABLE  4  (CONT. ) 


MACH  NUMBER  =  10  MACH  NUMBER  =  15  MACH  NUMBER  =  20 


DISTANCE 
FROM 
THROAT 
(  FT  ) 

0.05 

0.06 

0.07 

0.08 

0.09 

0.10 

0.11 

0.12 

0.13 

0.14 

0.15 

0.16 

0.17 

0.18 

0,19 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

0.80 

0.90 

1.00 

1.10 

1.20 

1.30 

1.40 

1.50 
1.60 

1.70 
1.80 
1.90 
2.00 
2.10 
2.20 

2.30 

2.40 

2.50 
2.60 

2.70 

2.80 


RADIUS 
OF  WALL 
(FT) 

0.10035 
0.10043 
0.10052 
0.10062 
0.10073 
0.10085 
0. 10098 
0.10112 
0.10127 
0.10143 
0.10160 
0.10178 
0.10197 
0.10216 
0.10237 
0.10258 
0.10516 
0.10847 
0.11246 
0.11705 
0.12221 
0 . 1 2  7°  0 
0.13409 
0.14077 
0.14791 
0.15549 
0.16352 
0.17199 
0.18087 
0.19015 
0.19985 
0.20994 
0.22041 
0.23127 
0.24251 
0.25411 
0.26607 
0.27841 
0.29107 
0.30408 
0.31740 
0.33105 


RADIUS 
OF  CORE 
(FT  ) 

0.10011 

0.10018 

0.10026 

0.10035 

0.10046 

0.10057 

0.10069 

0.10082 

0.10096 

0.10111 

0.10127 

0.10144 

0.10162 

0.10181 

0.10200 

0.10221 

0.10469 

0.10791 

0.11179 

0.11626 

0.12129 

0.12683 

0.13286 

0.13935 

0.14630 

0.15368 

0.16147 

0.16967 

0.17827 

0.18726 

0.19663 

0.20638 

0.21649 

0.22696 

0.23779 

0.24896 

0.26046 

0.27231 

0.28446 

0.29694 

0.30969 

0.32276 


RADIUS 
OF  WALL 
(  FT) 

0.03925 
0.04000 
0.04085 
0.04177 
0.04278 
0.04387 
0.04503 
0.04626 
0.04756 
0.04893 
0.05037 
0.05186 
0.05342 
0.05503 
0.05669 
0.05839 
0.07699 
0.09677 
0. 1] 688 
0.13709 
0.15735 
0 .17762 
0 .19792 
0.21823 
0.23859 
0.25897 
0.27941 
0.29988 
0.32039 
0.34066 
0.36064 
0.38033 
0.39973 
0.41882 
0.43764 
0.45617 
0.47442 
0.49242 
0.51014 
0.52762 
0.54483 
0.56182 


RADIUS 
OF  CORE 
(  FT  ) 

0.03915 
0.03989 
0 .04072 
0.04163 
0.04263 
0.04370 
0 . 04484 
0.04606 
0,04734 
0.04869 
0.05010 
0.05158 
0.05311 
0.05469 
0.05632 
0.05800 
0.07626 
0.09559 
0.11516 
0.13473 
0.15425 
0.17373 
0.19315 
0.21255 
0.23192 
0.25128 
0.27064 
0.28999 
0.30935 
0.32861 
0.34760 
0.36632 
0.38475 
0.40288 
0.42073 
0.43830 
0.45559 
0.47261 
0.48936 
0.50587 
0.52211 
0.53812 


RADIUS 
OF  WALL 
(  FT  ) 

0.01716 
0.01842 
0.01985 
0.02143 
0.02317 
0.02503 
0.02699 
0.02901 
0.03108 
0.03318 
0.03530 
0.03744 
0.03959 
0.04174 
0.04390 
0.04606 
0  a  Obi  55 
0.08882 
0.10996 
0.13109 
0.15226 
0.17349 
0.19471 
0.21538 
0.23559 
0.25535 
0.27467 
0.29357 
0.31209 
0.33024 
0.34803 
0.36550 
0.38264 
0.39949 
0.41604 
0,43232 
0.44834 
0.46410 
0.47963 
0.49492 
0.51000 
0.52485 


RADIUS 
OF  CORE 
(  FT  ) 

0.01707 
0.01831 
0.01971 
0.02127 
0.02297 
0.02479 
0.02671 
0.02869 
0.03072 
0.03277 
0.03484 
0.03692 
0.03901 
0.04110 
0.04319 
0.04528 
0.06595 
0.08615 
0.10602 
0.12571 
0.14529 
0.16482 
0.18429 
0.20345 
0.22215 
0.24039 
0.25820 
0.27560 
0.29259 
0.30921 
0.32548 
0.34140 
0.35701 
0.37231 
0.38731 
0.40203 
0.41649 
0.43068 
0 .44464 
0.45834 
0.47183 
0.48510 
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TABLE  4  ( CONT. ) 


MACH  NUMBER  =  10  MACH  NUMBER  =  15  MACH  NUMBER  =  20 


DISTANCE 
FROM 
THROAT 
(  FT) 

2.90 
3.00 

3.10 

3.20 

3.30 

3.40 

3.50 

3.60 

3.70 

3.80 

3.90 
4.00 

4.10 

4.20 

4.30 

4.40 

4.50 

4.60 

4.70 

4.80 

4.90 
5.00 

5.20 

5.40 

5.60 

5.80 

6.00 

6.20 

6.40 

6.60 

6.80 
7.00 

7.20 

7.40 

7.60 

7.80 

8.00 

8.20 

8.40 

8.60 

8.80 
9.00 


RADIUS 
OF  WALL 
(FT) 

0.34497 
0.35921 
0.37370 
0.38846 
0.40345 
0.41865 
0.43407 
0.44966 
0.46539 
0  ®  48 1  29 
0.49730 
0.51339 
0.52958 
0.54583 
0.56211 
0.57842 
0.59476 
0.61108 
0.62738 
0.64366 
0.65991 
0.67609 
0.70829 
0.74020 
0.77175 
0.80296 
0.83370 
0.86408 
0.89398 
0.92342 
0.95240 
0.98088 
1.00890 
1.03643 
1.06348 
1.09008 
1 .11619 
1.14186 
1 .16710 
1.19188 
1.21625 
1.24018 


RADIUS 
OF  CORE 
(FT  ) 

0.33607 
0.34968 
0.36351 
0.37758 
0.39188 
0.40635 
0.42104 
0.43587 
0.45083 
0 .46595 
0.48117 
0  ®  49646 
0.51184 
0.52728 
0.54274 
0.55823 
0.57375 
0.58924 
0.60471 
0.62017 
0.63559 
0.65095 
0.63150 
0.71177 
0.74170 
0.77128 
0.80043 
0.82919 
0.85748 
0.88533 
0.91272 
0.93961 
0.96605 
0.99200 
1.01747 
1.04248 
1.06701 
1.09109 
1 .11473 
1 .13790 
1.16065 
1.18297 


RADIUS 
OF  WALL 
(  FT  ) 

0.57857 
0 .59509 
0.61139 
0.62748 
0.64335 
0.65902 
0.67450 
0.68978 
0 .70487 
0.71979 
0.73452 
0.74908 
0.76346 
0.77769 
0.79175 
0.80564 
0.81939 
0 .83299 
0.84644 
0.85974 
0.87289 
0.88592 
0.91156 
0.93668 
0.96130 
0.98542 
1.00910 
1 .03234 
1.05514 
1 .07752 
1.09948 
1.12104 
1.14225 
1 .16309 
1 .18358 
1.20372 
1.22353 
1.24303 
1.26220 
1.28106 
1 .29964 
1.31791 


RADIUS 
OF  CORE 
(FT) 

0.55390 
0.56943 
0.58475 
0.59986 
0.61474 
0.62943 
0.64392 
0.65820 
0.67229 
0.68621 
0.69995 
0.71350 
0.72687 
0.74009 
0.75314 
0.76602 
0.77875 
0.79133 
0.80375 
0.81602 
0.82815 
0.84014 
0.86371 
0.88675 
0.90927 
0.93129 
0.95285 
0.97393 
0.99458 
1.01480 
1.03460 
1.05401 
1.07302 
1.09165 
1.10992 
1.12783 
1 .14539 
1.16263 
1.17953 
1 .19611 
1.21239 
1.22836 


RADIUS 
OF  WALL 
(FT) 

0.53949 
0*55394 
0.56819 
0.58225 
0.59613 
0.60984 
0.62338 
0.63675 
0.64997 
0.66302 
0.67593 
0.68870 
0.70131 
0.71379 
0.72614 
0.73836 
0.75044 
0.76240 
0.77424 
0.78596 
0.79757 
0 • 80905 
0.83170 
0.85392 
0.87573 
0.89717 
0.91821 
0.93890 
0.95925 
0.97925 
0.99894 
1.01831 
1.03739 
1  .056  18 
1.07469 
1.09293 
1 .11091 
1 .12863 
1 . 14610 
1 .16334 
1.18037 
1.19714 


RADIUS 
OF  CORE 
(FT) 

0.49815 

0.51100 

0.52365 

0.53611 

0.54839 

0.56049 

0.57241 

0.58417 

0.59576 

0.60719 

0.61847 

0.62961 

0.64059 

0.65143 

0.66213 

0.67271 

0.68314 

0.69344 

0.70363 

0.71369 

0.72363 

0.73345 

0.75275 

0.77162 

0.79006 

0.80810 

0.82574 

0.84302 

0.85994 

0.87650 

0.89274 

0*90864 

0.92423 

0.93952 

0.95451 

0.96921 

0.98365 

0.99780 

1.01170 

1.02534 

1.03873 

1.05188 
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TABLE  4  ( CONT. ) 


MACH  NUMBER  =  10  MACH  NUMBER  =  15  MACH  NUMBER  =  20 


DISTANCE 

FROM 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

RADIUS 

THROAT 

OF  WALL 

OF  CORE 

OF  WALL 

OF  CORE 

OF  WALL 

OF  CORE 

(  FT) 

(  FT  ) 

(FT  ) 

(  FT  ) 

(FT) 

(  FT  ) 

(  FT  ) 

9.20 

1.26369 

1.20486 

1.33591 

1.24404 

1.21370 

1.06481 

9.40 

1.28680 

1.22634 

1.35363 

1.25943 

1.23003 

1.07749 

9.60 

1.30953 

1.24743 

1.37108 

1.27454 

1.24612 

1.08996 

9.80 

1.33184 

1.26810 

1.38827 

1.28937 

1.26204 

1.10221 

10.0 

1.35378 

1.28839 

1.40520 

1.30394 

1.27777 

1.11425 

11.0 

1.45797 

1 .38422 

1.48624 

1.37298 

1.35359 

1.17141 

12.0 

1.55363 

1.47135 

1.56172 

1.43618 

1.42510 

1.22395 

13.0 

1.64153 

1.55056 

1,63216 

1.49409 

1.49273 

1.27231 

14.0 

1.72237 

1.62255 

1.69805 

TV  5472  2 

1.55684 

1.31690 

15.0 

1.79674 

1 .68795 

1 .75977 

1.59597 

1.61776 

1.35805 

16.0 

1.86518 

1.74730 

1.81766 

1.64070 

1.67574 

1.39604 

17.0 

1.92817 

1.80109 

1.87200 

1.68174 

1.73098 

1.43112 

18.0 

1.98614 

1.84978 

1.92307 

1.71937 

1.78369 

1.46353 

19.0 

2.03944 

1.89373 

1.97108 

1.75382 

1.83402 

1.49344 

20.0 

2.08844 

1.93334 

2.01621 

1.78532 

1.88212 

1.52103 

21.0 

2.13341 

1.96891 

2.05865 

1.81408 

1.92812 

1.54646 

2  2.0 

2.17465 

2.00074 

2.09857 

1.84027 

1.97212 

1.56987 

23.0 

2.21242 

2.02Q12 

2.13608 

1.86406 

2.01422 

1.59139 

24.0 

2.24693 

2  .  n  5  4  3  0 

2.17134 

1  .88560 

2.05451 

1.61113 

2  5.0 

2.27842 

2 .07652 

2.20446 

1.90504 

2.09307 

1.62920 

26.0 

2.30708 

2.09600 

2.23554 

1.92250 

2.12997 

1.64571 

27.0 

2.33309 

2.11295 

2.26469 

1.93812 

2.16528 

1.66074 

28.0 

2.35663 

2.12757 

2.29199 

1.95199 

2.19905 

1.67437 

29.0 

2.37787 

2 .14006 

2.31754 

1.96424 

2.23135 

1.68670 

30.0 

2.39696 

2.15059 

2.34140 

1.97496 

2.26222 

1.69778 

31.0 

2.41405 

2.15935 

2.36366 

1.98425 

2.29170 

1.70771 

32.0 

2.42928 

2.16650 

2.38437 

1.99220 

2.31984 

1.71653 

33.0 

2.44279 

2.17221 

2.40360 

1.99891 

2.34666 

1.72432 

34.0 

2.45471 

2.17664 

2.42141 

2 . 00446 

2.37222 

1.73113 

35.0 

2.46516 

2.17996 

2.43784 

2.00893 

2.39652 

1.73703 

36.0 

2.47428 

2.18233 

2.45293 

2.01242 

2.41960 

1.74207 

37.0 

2.48220 

2.18391 

2.46672 

2.01501 

2.44148 

1.74631 

3  8.0 

2.48903 

2 .18486 

2.47921 

2.01679 

2.46217 

1.74980 

38.2 

2.49028 

2.18499 

2.48155 

2.01706 

2.46617 

1.75041 

38.4 

2.49149 

2.18510 

2.48384 

2.01730 

2.47012 

1.75100 

38.6 

2.49266 

2.18520 

2.48607 

2.01751 

2.47402 

1.75155 

38.8 

2.49380 

2.18528 

2.48825 

2.01770 

2.47787 

1.75209 

39.0 

2.49491 

2.18535 

2.49038 

2.01786 

2.48168 

1.75259 

39.2 

2.49599 

2.18541 

2.49244 

2.01800 

2.48544 

1.75307 

39.4 

2.49703 

2.18545 

2.49445 

2.01812 

2.48915 

1.75352 

39.6 

2.49805 

2.18549 

2.49639 

2.01821 

2.49282 

1.75396 

39.8 

2.49904 

2.18552 

2.49825 

2.01828 

2.49643 

1.75436 

o 

9 

o 

2.50000 

2.18554 

2 .50000 

2.01834 

2 .50000 

1.75474 
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TABLE  5 


CARD  INPUT  DATA  FOR  ACP*  AND  RGTBLP**  PROGRAMS 


MACH 

NUMBER 

10 

-  ACP* 

NN 

IFTH 

LM 

LN  KK 

EMT 

XT  RCT 

EMZP 

10 

0 

2 

0  100 

10. 

9  .  0  . 

2.8 

GAMA 

XZ 

DZX 

DTX 

EPSA 

EPSB 

1 

. 4E+00 

1 .0E-03 

2. 0E-03 

2.0E-03 

5.0E-06 

5.0E-06 

MACH 

NUMBER 

10 

-  RGTBLP** 

THETA 

REXIT  XPUNCH 

XO 

RO 

OMEGAO 

HI 

.000 

3  2 

.185577  100. 

-.75 

.33924 

-5. 

0. 

NN 

LM 

NFS 

IP  J 

LAWCF 

910 

2 

1 

1  2 

1 

M 

TX  (  2  ) 

TX ( 3 )  TX ( 4 ) 

TX ( 5 )  TX ( 6 ) 

TX  (  8  ) 

2 

900. 

2.0 

700. 

10.  400. 

300. 

DELX 

XEND 

IP 

.01 

- 

.28 

.001 

.02 

10 

.01 

0.2 

.05 

5.0 

2 

.2 

999. 

MACH 

NUMBER 

15 

-  ACP* 

NN 

IFTH 

LM 

LN  KK 

EMT 

XT  RCT 

EMZP 

15 

2 

0  100 

15. 

0.  0. 

0  © 

GAMA 

XZ 

DZX 

DTX 

EPSA 

EPSB 

1 

•  4E  +  00 

1 .0E-03 

2 .0E-03 

2.0E-03 

5.0E-06 

5  •  OE— 06 

OMEGA  EMB  XER  RCPLUS 

11.0  4.134  0.99  12.0 


*  AXISYMMETRIC  CORE  PROGRAM  (SEE  REFERENCE  (15)) 

**  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM  (SEE  APPENDIX  C) 
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TABLE  5  (CONT.) 

MACH  NUMBER  15  -  RGTSLP** 

THETA  REXIT  XPUNCH  XO  RO  OMEGAO  HI 

.0003  2.018408  100.  -.45  .1914  -3.  0. 

NN  LM  NFS  IP  J  LAWCF 

<515  2  112  1 

M  TX ( 2 )  T  X ( 3 )  T  X ( 4  >  T  X ( 5 )  TX(6)  TX(8) 

2  1200.  2.0  900.  10.  400.  300. 


DELX 

XEND 

IP 

.01 

-.28 

.001 

.02 

10 

.01 

0.2 

.05 

5.0 

2 

.2 

999  . 

MACH  NUMBER  20  -  ACP* 


I  FTH 

LM 

LN 

KK 

EMT 

XT 

RCT 

EMZP 

2 

2 

0 

100 

20. 

0. 

0  e 

0. 

GAMA 

xz 

DZX 

DTX 

EPSA 

1.4E+00 

1  ©  0  E“" 

■03 

2 • OE-O  3 

2.0E-03 

5 . 0E-06 

OMEGA  EMB  XER  RCPLUS 

11.0  4.534  0.99  16.0 


MACH  NUMBER  20  -  RGTBLP** 

THETA  REXIT  XPUNCH  XO  RO  OMEGAO  HI 

.0001  1.758123  100.  -.4  .16657  -5.  0. 

NN  LM  NFS  IP  J  LAWCF 

Q  2  0  2  1  1  2  1 

M  TX ( 2 )  T X ( 3 )  T X ( 4 }  T X ( 5 )  T X ( 6  )  TX ( 8  ) 

2  1600.  2.0  1200.  5.0  800.  300. 


DELX 

XEND 

IP 

.01 

-.28 

.001 

.02 

10 

.01 

0.2 

.05 

5.0 

2 

.2 

999. 

*  AXISYMMETRIC  CORE  PROGRAM  (SEE  REFERENCE  (15)) 

**  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM  (SEE  APPENDIX  C) 


EPSB 

5.0E-06 


! 
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TABLE  6 

SELECTED  NOZZLE  DIMENSIONS 


MACH  NUMBFR 

1  0 

1  5 

20 

LENGTH  (FT! 

) 

INLET 

0.75 

0.45 

0.40 

EXPANSION 

40.00 

40.00 

40.00 

TOTAL 

40.75 

40.45 

40.40 

DIAMETER  (FT) 

ENTRANCE 

0.679 

0.383 

0.333 

THROAT 

0.200 

0.074 

0.028 

EXIT 

5.000 

5.000 

5.000 

CORE* 

4.371 

4.037 

3.509 

USABLE** 

4.005 

3.752 

3.282 

AREA  (FT2) 

ENTRANCE 

0.3622 

0.1154 

0.0872 

THROAT 

0.0315 

0.0043 

0.0006 

EXIT 

19.635 

19.635 

19.635 

CORE* 

15.006 

12.798 

9.6729 

USABLE** 

12.596 

11.059 

8.4588 

*  INVISCID 

CORE  MEASURED  AT  NOZZLE  EXI' 

**  USABLE  DTAMFTFR  AT 

NOZZLE  EXIT 

DEFINI 

AS  THE  NOZZLE  EXIT  DIAMETER  MINUS  TWICE 
THE  BOUNDARY  LAYER  THICKNESS  AT  THE  EXIT 
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TABLE  7 

CORRECTED  WALL  RADII  IN  SUBSONIC  INLETS 


MACH 

NUMBER 

POINT 

NUMBER 

DISTANCE 

FROM 

THROAT 

(FEET) 

DISTANCE 
FROM 
THROAT 
( INCHES) 

RADIUS 

OF  WALL 
(INCHES) 

RADIUS 
OF  WALL 
( FEET ) 

10 

37 

-0.39 

-4.68 

2.166 

0. 18050 

38 

-0.38 

-4.56 

2.108 

0.17567 

39 

-0.37 

-4.44 

2.052 

0.17100 

40 

-0.36 

-4.32 

1.996 

0.16633 

15 

58 

-0.24 

-2.88 

1.259 

0.10492 

68 

-0.23 

-2.76 

1.196 

0.09967 

78 

-0.22 

-2.64 

1.136 

0.09467 

88 

-0.21 

-2.52 

1.077 

0.08975 

98 

-0.20 

-2.40 

1.021 

0.08508 

108 

-0.19 

-2.28 

0.968 

0.08067 

118 

-0.18 

-2.16 

0.918 

0.07650 

33 

-0.26 

-3.12 

1.392 

0. 11600 

43 

-0.25 

-3.00 

1.327 

0.11058 

53 

-0.24 

-2.88 

1.262 

0.10517 

20 
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TEMPERATURE  (°K) 


FIG.  1  COMPARISON  OF  PRESSURE -TEMPERATURE  VARIATIONS  AT  HIGH  DENSITIES  WITH 
AEDC  DATA  TABULATIONS 


REYNOLDS  NUMBER  PER  FOOT 
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FIG.  2 


VARIATION  OF  REYNOLDS  NUMBER  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONSTANT  PRESSURE  OR  CONDENSATION  THRESHOLD  OPERATION  AT  MACH  20 
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FIG.  3 


VARIATION  OF  PRESSURE  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 


DENSITY  (AMAGATS) 


NOLTR  71-6 


40  44  48  52  56  60 

TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG.  4  VARIATION  OF  DENSITY  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 


ENTROPY  (S/R) 
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40  44  48  52  56  60 

TEMPERATURE  IN  TEST  SECTION  (°K) 

FIG.  5  VARIATION  OF  ENTROPY  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 
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40  44  48  52  56  60 

TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG.  6  VARIATION  OF  ENTHALPY  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 


SOUND  SPEED  (FT/SEC) 
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TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG  7  VARIATION  OF  SOUND  SPEED  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 


MASS  FLOW  PER  UNIT  AREA/MACH  NUMBER  (LBM/FT2  SEC 
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FIG.  8 


40  44  48  52  56  60 

TEMPERATURE  IN  TEST  SECTION  (°K) 


VARIATION  OF  MASS  FLOW  IN  TEST  SECTION  WITH  TEMPERATURE  FOR 
CONDENSATION  THRESHOLD  OPERATION 
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TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG.  9  VARIATION  OF  REYNOLDS  NUMBER  IN  TEST  SECTION  WITH  TEMPERATURE 
FOR  CONDENSATION  THRESHOLD  OPERATION 


VISCOSITY  (LBAVFTSEC) 
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40  44  48  52  56  60 

TEMPERATURE  IN  TEST  SECTION  (°K) 

FIG.  10  VARIATION  OF  VISCOSITY  IN  TEST  SECTION  WITH  TEMPERATURE 
FOR  CONDENSATION  THRESHOLD  OPERATION 


REDUCED  VISCOSITY  AT  ATMOSPHERIC  PRESSURE,  /**//*’ 
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in 


FIG.  11  VISCOSITY-TEMPERATURE  VARIATION  USED  IN  BOUNDARY 
LAYER  CALCULATIONS 


SUPPLY  PRESSURE  (ATM) 
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40  44  48  52  56  60 


TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG.  12  VARIATION  OF  SUPPLY  PRESSURE  REQUIRED  FOR  CONDENSATION  THRESHOLD 
OPERATION  WITH  TEST  SECTION  TEMPERATURE  AND  MACH  NUMBER 


SUPPLY  TEMPERATURE  (°K) 
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FIG.  13  VARIATION  OF  SUPPLY  TEMPERATURE  REQUIRED  FOR  CONDENSATION  THRESHOLD 
OPERATION  WITH  TEST  SECTION  TEMPERATURE  AND  MACH  NUMBER 


SUPPLY  DENSITY  (AMAGATS) 
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TEMPERATURE  IN  TEST  SECTION  (°K) 


FIG.  14  VARIATION  OF  SUPPLY  DENSITY  REQUIRED  FOR  CONDENSATION  THRESHOLD 
OPERATION  WITH  TEST  SECTION  TEMPERATURE  AND  MACH  NUMBER 


SUPPLY  PRESSURE  (ATM) 
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FIG.  15  VARIATION  OF  SUPPLY  PRESSURE  REQUIRED  FOR  CONDENSATION  THRESHOLD 

OPERATION  WITH  TEST  SECTION  REYNOLDS  NUMBER  PER  FOOT  AND  MACH  NUMBER 
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RE 


FIG.  16  REQUIRED  SUPPLY  PRESSURE,  CORE  DIAMETER,  AND  VOLUME  FLOW  RATE  FOR  GIVEN 
REDc  AND  MASS  FLOW  RATE  (MACH  10) 
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m  (LBM/SEC) 


FIG.  17  REQUIRED  SUPPLY  PRESSURE,  CORE  DIAMETER,  AND  VOLUME  FLOW  RATE  FOR  GIVEN 
REDc  AND  MASS  FLOW  RATE  (MACH  15) 
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m  (LBM/SEC) 


FIG.  18  REQUIRED  SUPPLY  PRESSURE,  CORE  DIAMETER,  AND  VOLUME  FLOW  RATE  FOR  GIVEN 
REDc  AND  MASS  FLOW  RATE  (MACH  20) 


NOZZLE  WALL  TEMPERATURE  (°K) 
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FIG.  19  ASSUMED  VARIATION  OF  NOZZLE  WALL  TEMPERATURE  WITH  AXIAL  DISTANCE 
DOWNSTREAM  OF  THROAT 


BOUNDARY  LAYER  MOMENTUM  THICKNESS 
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CUTOFF 

POINT 


USABLE 
>  TEST 
RADIUS 


FIG.  21  SCHEMATIC  DIAGRAM  OF  NOZZLE  EXIT  SHOWING  DEFINITION  AND  LOCATION  OF 
OPTIMUM  CUTOFF  POINT 


FIG.  22  VARIATION  OF  SURFACE  HEAT  TRANSFER  RATE  WITH  AXIAL  DISTANCE  IN 
THROAT  REGION  FOR  ASSUMED  TEMPERATURE  DISTRIBUTIONS 


FIG.  23  VARIATION  OF  CALCULATED  ANGLE  AT  INFLECTION  POINT  WITH  INPUT  VALUE 
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CALCULATED  ANGLE  AT  INFLECTION  POINT  (DEGREE) 


FIG.  24  VARIATION  OF  CALCULATED  ANGLE  AT  INFLECTION  POINT  WITH  THROAT 
RADIUS  OF  CURVATURE/THROAT  RADIUS 
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APPENDIX  A 

MODIFIED  CENTERLINE  MACK  NUMBER  DISTRIBUTION 


The  Axisyraraetric  Core  Program  CACP)  has  been  used  to  compute 
the  inviscid  core  flow  in  the  supersonic  portion  of  the  nozzle.  This 
program  has  been  modified  from  the  form  presented  in  reference  C15) 
in  order  to  incorporate  a  third  centerline  Mach  number  distribution. 
For  high  Mach  numbers  where  throat  diameters  are  typically  small ,  the 
cubic  and  exponential  distributions  contained  in  the  ACP  program  lead 
to  large  radii  of  curvature  at  the  nozzle  throat.  Since  such  nozzles 
are  difficult  to  manufacture  and  can  lead  to  throat  heating  problems, 
a  third  centerline  Mach  number  distribution  has  been  added  to  the 
program. 

In  this  third  distribution,  a  quadratic  Mach  number  variation  in 
the  region  of  the  nozzle  throat  is  matched  to  one  corresponding  to  a 
region  of  radial  flow  further  downstream.*  However,  the  Mach  number 
in  this  radial  flow  region  does  not  reach  the  exit  Mach  number  with 
the  desired  zero  gradient.  To  smooth  this  transition,  a  section  of 
variable  length  is  inserted  between  the  radial  flow  region  and  the 
uniform  flow  region.  In  this  section,  the  Mach  number  is  given  by 
the  expression 


..  . .  “  a  ( x ,  “  x )  .  i 

M  =  Mg  e  t  (A-l) 

where  Mg  is  the  exit  Mach  number  and  a  and  x^  are  determined  by 
matching  the  Mach  number  and  its  gradient  at  some  arbitrary  point  to 
the  values  obtained  from  the  radial  flow  region  distribution.  This 
match  point  then  becomes  the  boundary  between  the  radial  flow  region 
and  the  transition  region.  The  match  point  is  chosen  by  specifying 
*See  the  derivation  of  the  Mach  number  distribution  on  pages  A- 3  to  A-7. 
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its  location  as  some  fraction  (XER)  of  the  distance  from  the  throat 
at  which  the  exit  Mach  number  would  be  reached  if  no  transition 
section  were  inserted.  The  quantity  xt  is  the  location  at  which  the 
transition  section  ends  and  the  uniform  flow  region  begins. 

The  use  of  this  mpltiple-section  distribution  will  decrease  the 
radius  of  curvature  at  the  nozzle  throat  while  allowing  the  Mach 
number  to  vary  smoothly  along  the  nozzle  axis.  It  has  been  incorporated 
into  the  program  and  identified  as  the  Harris  centerline  Mach  number 
distribution  to  distinguish  it  from  the  Schaaf  (exponential)  and 
Thicks tun  (cubic)  distributions.  The  Harris  centerline  is  used  if  the 
input  quantity  IFTH  equals  2.  In  order  that  the  instructions  given  in 
reference  (15)  for  running  this  program  with  the  Schaaf  and  Thickstun 
centerline  distributions  remain  valid,  the  input  quantities  required 
for  defining  the  Harris  distribution  are  read  from  a  third  data  card 
(when  IFTH  equals  2)  with  the  statements 

45  READ  (5,  46)  OMEGA,  EMB,  XER,  RCPLUS 

46  FORMAT  (4F10.5) 

The  required  input  quantities  are  defined  as  follows  (see  sketch  below) : 

OMEGA  angle  at  the  inflection  point  measured  with  respect 
to  the  nozzle  centerline  (see  next  section) , 

EMB  Mach  number  at  point  B, 

XER  location  of  the  match  point  between  the  radial  flow 

region  and  the  transition  section  expressed  as  a 
fraction  of  the  distance  from  the  throat  at  which  the 
exit  Mach  number  would  be  reached  if  no  transition 
section  were  inserted, 

RCPLUS  desired  ratio  of  throat  radius  of  curvature  to  throat 
radius 

The  values  of  XT,  RCT,  and  EMZP  (defined  in  reference  (15))  are 
calculated  in  subroutine  CLP.  Therefore,  dummy  values  (or  blanks)  can 
be  used  for  these  three  quantities  on  the  first  data  card.  The  values 
of  XT,  RCT,  and  EMZP  written  in  the  output  heading  are  the  calculated 
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ones.  The  input  quantities  OMEGA,  EMB ,  and  MER  a, ye  written  in  the 
output  as  PA,  PB,  and  PC,  respectively,  while  the  coefficient  a  in  the 
transition  region  Mach  number  distribution  (equation  A-l)  is  written 
as  PD.  The  value  input  for  the  quantity  RCPLUS  is  not  written  in  the 
output  but  can  be  obtained  by  dividing  the  output  values  of  RCT  and 
THROAT  RADIUS. 

DERIVATION  OF  MACH  NUMBER  DISTRIBUTION 


The  derivation  which  follows  is  strictly  valrd  only  when  the  ratio 
of  specific  heats  y  is  constant  with  a  value  of  1.4  throughout  the 
flow.  Thus,  for  flows  with  other  constant  or  variable  specific  heats, 
the  flow  will  not  be  completely  radial  in  the  region  BCD  shown  in  the 
above  sketch.  However,  this  will  not  cause  any  real  difficulty, 
because  the  flow  in  the  region  BCD  is  calculated  from  the  resulting 
centerline  Mach  number  distribution  using  the  method  of  characteristics 
rather  than  from  the  radial  flow  equations.  In  other  words,  the  radial 
flow  region  BCD  exists  formally  only  in  the  following  derivation.  It 
does  not  exist  formally  in  the  nozzle  calculation.  As  the  gas  departs 
from  perfect  gas  behavior,  not  only  does  the  flow  in  region  BCD  depart 
somewhat  from  strictly  radial  flow,  but  also  the  calculated  angle  at 
the  inflection  point  departs  from  the  input  quantity  OMEGA  used  to 
determine  the  centerline  Mach  number  distribution.  Figure  23  shows 
that  for  a  given  RCPLUS  (ratio  of  throat  radius  of  curvature  to  throat 
radius) ,  the  calculated  inflection  point  angle  matches  the  input  value 
OMEGA  fairly  well  at  larger  values  of  OMEGA,  but  departs  more  and  more 
as  OMEGA  decreases.  Moreover,  for  the  case  shown  in  Figure  23,  the 
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minimum  value  obtained  for  the  calculated  angle  is  aboye  the  desirable 
10°-12°  range.  Figure  24  shows  that  by  increasing  RCPLTJS  the 
calculated  angle  can  be  made  to  agree  more  closely  with  the  input 
value  OMEGA  and  (as  was  also  shown  in  Figure  23)  the  agreement  is 
better  at  the  larger  values  of  OMEGA.  However,  it  is  not  necessary  to 
make  the  two  angles  agree  closely.  The  important  thing  is  that  the 
calculated  angle  be  kept  within  the  desired  range.  As  the  two  figures 
have  shown,  this  may  require  the  use  of  a  value  of  OMEGA  substantially 
smaller  than  the  desired  angle  at  the  inflection  point.  At  the  present 
time,  these  figures  are  about  the  only  guide  to  selecting  the  value  of 
OMEGA.  Until  more  such  results  are  available,  the  process  of  obtaining 
the  desired  angle  will  be  a  trial  and  error  one. 

In  the  following  derivation,  the  sonic  line  is  assumed  to  be 
straight,  perpendicular  to  the  nozzle  axis,  and  located  at  point  A. 

The  Mach  number  distribution  from  A  to  B  is  quadratic  in  distance  with 
the  same  slope  at  B  as  that  calculated  for  point  B  using  the  radial 
flow  assumptions  for  region  BCD.  The  following  notation  is  used. 


r 

e 

x  (=x  r  ) 
e 

y  (=y  re) 

M 


D 


M  ' 

B 


radius  of  inviscid  core  at  exit 

distance  from  throat 

distance  from  apparent  source  of  radial  flow 
(point  0) 

Mach  number  at  point  D 

ciM 

Mach  number  gradient  — —  at  point  B 

dx 


From  A  to  B 


M(x+) 


MB  “  MA  „+ 
2  XB 


«)’*■;  -;(S) 


M, 


(A- 2) 


A- 4 
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-  M^) 


+ 


+  M 


B 


A 


where 


it 


X 


B 


M 


A 


B 


mb  +  ma 


(A- 3) 


(A- 4) 


Since  point  A  is  the  position  of  the  sonic  line,  M  equals  1.  The 

i  A 

gradient  M  at  the  sonic  point  is  found  f^om  the+throat  radius  r* 

1  c  ^c 

and  the  input  quantity  RCPLUS  (= —  =  — qr)  using  the  perfect 


,  + 
(=r*r 


gas  relation 


i 


M 


M=1 


1_  fl.2 

+  \  RCPLUS 
r*  * 


(A- 5) 


From  B  to  D 

For  the  radial  flow  region  BCD,  the  flow  angle  to  at  the  inflection 
point  c  is  related  to  the  Mach  numbers  at  points  B  and  D.  The  relation¬ 
ship  is 


v  (Mg)  +  4to  =  v  (Mg) 


(A- 6) 


where  v (M)  is  the  Prandtl-Meyer  angle  tabulated  for  two-dimensional 

flow,  (For  table  of  v  (M) ,  see  reference  (27)).  Thus,  the  Mach  number 

at  point  B  (input  quantity  EMB)  is  obtained  from  the  Mach  number  at 

point  D  (=  design  Mach  number)  and  the  desired  maximum  expansion  angle 

to  (input  quantity  OMEGA).  Typically,  to  is  in  the  range  10°  -  12°. 

In  this  radial  flow  region,  the  distance  y  from  the  source 

(point  0)  is  related  to  the  Mach  number  by  the  one-dimensional  relation 

for  the  area  ratio  —  , 

A* 
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1 

M 


.  2 
■  Y+l 


Cl  +  M2}]  2tY  1} 


where  y*  is  the  value  of  y  for  which  M  =  1.  For  y  =  1.4,  this 
written  as 


=  IK  CM)  ]  1/2 

y  * 


where 


Also 


where 


Noting  that 


K  CM) 


1  .1  +  0. 2M2  3 
M  L  1.2  J 


—  =  g  (u> ,  Md) 
e 


gt^Hp)  “  2(sin  ?>  tKWD))l/2 


X  =  y  -  OA  =  y  -  CyB  -  xB) 


we  find 


1 

x+  =  —  =  g(o),MD)  {  [K  (M)  ]  2 

e 


[k(mb)]2  }  +  p- 

e 


and  by  differentiation 


dM  _  2. 4M 

dx+  g(.03,MD)  CM2-1) 


1 . 2M 
1+0. 2M2 


1 

2 


(A- 7) 


can  be 


(A-8) 


(A-9 ) 


(A-10) 


(A-ll) 


(A-12 ) 


(A-13) 
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SOLUTION  PROCEDURE 

The  procedure  to  be  used  in  setting  up  this  centerline  Mach 

number  distribution  is  now  as  follows.  First,  the  nozzle  design  Mach 

number  M  and  the  desired  maximum  expansion  angle  w  are  chosen.  Then 

using  equation  (A-6)  and  tables  such  as  in  reference  (27) ,  the  Mach 

number  at  point  B  (=M0)  is  calculated.  The  three  quantities  Mp,  to, 

M_  and  the  desired  ratio  of  the  throat  radius  of  curvature  to  throat 
B  p 

Xv-C 

radius  —  are  read  into  the  program  as  the  input  quantities  EMT, 

OMEGA,  EMB,  and  RCPLUS ,  respectively.  In  the  program,  using  the 
throat  radius  r*+  obtained  from  the  area  ratio  -  Mach  number  relation- 

it 

ship  for  one-dimensional  flow,  equations  (A- 5) ,  (A- 4) ,  (A- 9)  ,  (A-ll) , 

I 

and  (A-13)  are  used  in  the  order  given  to  calculate  the  quantities  M  , 

I  ^ 

x+,  K(M_),  K  (M  )  ,  g(oo,l{L),  and  M  ,  respectively.  With  these  quantities, 

ID  ID  U  JJ  ID 

M  and  M'  can  be  determined  from 


0<x<xD 

Jd 

xT  <x<XER*x 

JD  L) 

XER*xD<xt 


equations  (A-2)  and  (A-3) 
equations  (A-12)  and  (A-13) 
equation  (A-l) 


Since  equation  (A-12)  can  not  be  solved  explicitly  for  M,  a  Newton- 
Raphson  iterative  procedure  is  used.  Beginning  with  M  at  the  previous 
station  as  a  first  guess,  the  value  of  M  is  corrected  in  successive 
iterations  until  two  consecutive  values  agree  to  within  0.002  percent. 


*Ecraation  (A- 9)  is  used  twice,  once  for  KCM-),  once  for  K(M_). 

-D  D 
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LISTING  OF  MODIFIED  SUBROUTINES 


SUBROUTINE  CLP  ( EMZP , P A , PB » PC , PD , EMT , XT , RC T , I FTH * RSTAR , EPSB . GAMA ) 
C  ACP  1/30/71  WJG 

COMMON  HSTAG,SSTAG,PSTAG»RHOSTG»ASTAR»TA*TB*TC*TD*TE»TF*TG 
IF  (IFTH-l)  1,1*45 
1  IF(EMZP)  5,5,10 

5  EMZP=SQRT  ( ( ( GAMA+1 . ) /2 . ) /RSTAR/RCT ) 

10  IF(RCT)  15,15,20 

15  RC T  = ( (GAMA+1. )  /2, ) / ( RSTAR*EMZP**2 ) 

20  IF(IFTH)  30,30,25 
C  THICKSTUN  CENTER  LINE 

25  PA= ( EMT-1 • ) /XT**2 

PB=2./XT-EMZP/ ( EMT-1.  ) 

PC=0.0 
PD=0 . 0 

IFIPB+l./XT)  30,30,75 
C  SCHAAF  CENTER  LINE 

30  IFTH=0 
PB  =  0. 1 

DO  35  L  =  1 »  2  5 

PT  =  PB*XT 

EBX  =  EXP  (PT*XT) 

EMX  =  EMZP/ ( 2 • *P  T ) 

PC  =  EMX* ( 1 . -EBX )  +EBX* ( EMT  -  1.) 

PCP  =  EMX*( EBX-1 . ) /PB  +  EBX*XT**2*(EMT-1,-EMX> 

HB  =  PC/PCP 
PB  =  PB  -  HB 
PD  =  L 

I F ( A8S  (HB)  -  EPSB)  40,40,35 
35  CONTINUE 

40  PA  =  EMZP/( 2,*XT*PB)  +  1.  -  EMT 
GO  TO  75 

C  HARRIS  CENTER  LINE 

45  READ  (5,46)  OMEG A , EMB , XER , RCPLUS 

46  FORMAT  (4F10.5) 

TE= ( ( l.+.2*EMT**2 >/1.2) **3/ EMT 
TD= ( <1.+.2*EMB**2)/1.2)**3/EMB 
TC=l./( 2. *SIN( OMEGA/ 114.5916 )*SQRT(TE) ) 

TB  =  2.4*EMB*SQRT ( 1 . 2*EMB / ( 1 .  +  . 2*EMB**2 ) ) / ( TC* ( EMB**2-1 .  ) ) 
RCT=RCPLUS*RSTAR 

EMZP=SORT  ( ( (GAMA+1. ) /2 .) /RSTAR/RCT ) 

AXB=(TB-EMZP) /2. 

XB=(EMB-1.)/(AXB+EMZP) 

TA  =  XB 
PA=OMEGA 
PB=EMB 
PC=XER 

XT=TC*< SORT (TE) -SORT (TD) >+TA 
X  E  =  XER*XT 

TF= ( (XE-TA) /TC+SQRT (TD) )**2 
EM=EMT-1«5* (XT-XE) 

1=0 

50  1=1+1 

EMSQ=EM**2 
FM= ( 5.  +  EMSQ  )  /6 • 

DM= (EM*FM-TF*EMSQ/FM**2 ) / ( FM-EMSQ) 

EM=EM+DM 

IF  (ABS(DM) /EM-2.0E-05)  70,70,55 
55  IF  (1-25)  50,50,60 
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60  WRITE  (6,65)  EM, DM 

65  FORMAT  (26H0  NO  CONVERGENCE  IN  CLP  1P2E15.6//) 

70  EME=EM 

EMPE=3.456*SQRT (TF) * ( EM/ ( 1 .+. 2*EM**2 ) ) **2 / ( TC* ( EM**2- 1 . ) ) 

ADX=0.5*EMPE/EME 

ADXSQ=ALOG ( EMT /EME ) 

DX=ADXSQ/ADX 
A=ADX/DX 
XT  =  DX+X  E 
PD  =  A 
TG  =  XE 
TE=AXB/XB 
TF=EMZP 
75  RETURN 
END 


SUBROUTINE  CL  ( EMT , XT , PA , PB , PC , PD , EM , EMP , X , I FTH ) 

C  ACP  1/30/71  WJG 

COMMON  HSTAG»SSTAG»PSTAG,RHOSTG»ASTAR»TA,TB,TC,TD,TE»TF»TG 
IF  (IFTH-1)  10,5,15 
C  THICKSTUN  CENTER  LINE 

5  EM=EMT-PA*< l.+PB*X)*( X-XT >**2 
EMP  =  -PA*(X-XT)*(PB*(3.*X-XT  )+2.  ) 

GO  TO  60 

C  SCHAAF  CENTER  LINE 

10  EBX  =  EXP  (PB*(X-XT)**2) 

EM  =  EMT  -  PA*(EBX-1.) 

EMP  =  -2.*PA*PB*(X-XT)*EBX 
GO  TO  60 

C  HARRIS  CENTER  LINE 

15  1=0 

EMZP=TF 

IF  (X-TA)  20,20,25 
20  EM=1.+X*(EMZP+TE*X) 

EMP=EMZP+2.*TE*X 
GO  TO  55 

25  IF  (X-TG)  30,50,50 
30  FKM=( (X-TA) /TC+SQRTITD) )**2 
EM=EM1 
35  1=1+1 

EMSQ=EM**2 
FM= ( 5.  +  EMSQ  )  /6. 

DM= ( EM*FM-FKM*EMSQ/FM**2 ) / ( FM-EMSQ ) 

EM=EM+DM 

IF  ( ABS(DM) /EM-2.0E-05)  45,45,40 
40  IF  (1-25)  35,35,45 

45  EMP  =  3.456*SQRT (FKM)*(EM/ ( 1 .  +  . 2*EM**2 ) ) **2 / ( TC*  (  EM**2~1 .  ) ) 
GO  TO  55 
50  DX=XT-X 

EM=EMT*EXP ( -PD*DX**2 ) 

EMP  =  2.-!<'PD*DX*EM 
55  EM1=EM 
60  RETURN 
END 
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APPENDIX  D 

UTILITY  PROGRAM  FOR  PREPARING 
BOUNDARY  LAYER  INPUT  DATA  CARDS 


In  order  to  calculate  the  subsonic  inlet  contour  using  the  option 
available  in  the  boundary  layer  program,  the  gas  properties  for  Mach 
numbers  less  than  one  must  be  supplied  to  the  program.  These  properties 
are  read  into  the  boundary  layer  program  from  cards  to  avoid  modifying 
the  tape  from  the  Axisymmetric  Core  Program  which  contains  the  gas 
properties  only  for  the  supersonic  Mach  numbers.  If  the  Mach  number 
at  the  beginning  of  the  inlet  is  not  too  low,  then  the  gas  properties 
calculated  during  the  isentropic  expansion  calculation  may  be  suf¬ 
ficient  to  give  the  desired  accuracy  in  the  required  interpolations. 

If  not,  the  isentropic  expansion  can  be  recalculated  and  more  subsonic 
data  generated.  In  either  case,  the  subsonic  gas  property  data 
present  on  the  isentropic  expansion  tape  can  be  transferred  to  cards 
in  the  form  required  for  the  boundary  layer  program  by  using  the 
utility  program  SUBSON  listed  below. 

The  only  input  to  the  program  is  the  isentropic  expansion  tape 
which  is  mounted  on  unit  17.  The  program  will  read  a  set  of  sub¬ 
sonic  data  from  the  tape,  convert  it  to  the  correct  form,  and  punch 
the  results  into  cards,  repeating  this  procedure  until  the  area  ratio 
decreases  below  1.0001.  The  data,  for  area  ratio  equal  to  one  (throat 
or  sonic  conditions)  is  then  punched  into  the  final  card. 
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PROGRAM  LISTING 


1  READ ( 17*2)  NN.PSTAG.TT.RHOSTG.HSTAG.TZ.SSTAG.TAR.TSTAR.ASTAR. 
1PROSTR.RHOSTR*  RHOQST  *EMMAX»DELTAP 

2  FORMAT ( 37H1  ISENTROPIC  EXPANSION  TABLE  NUMBER  I4//18H  SUPPLY  COND 
1  IT  I ONS/9H  P(ATM.)=1PE14.7,11H  T < KEL V  I N ) =1PE1 4 . 7 , 22H  RHO ( SLUGS/CUB  I 
2C  FT. )  =  1PE14.7/23H  ENTHALPY ( BTU/LB . MASS ) =1PE14 . 7 * 4H  Z=1PE14.7, 

36H  S/R=1PE14«7»10H  A/ ASTAR= 1 PE14, 7/ / 1 8H  THROAT  COND I T I ONS/ 1 1H  T< 
4KELVIN) =1PE14.7»18H  ASTAR ( FT . /SEC . )=1PE14.7.10H  P(ATM.)=1PE14.7/ 

522H  RHO (SLUGS/CUBIC  F T .  )  =  1 PEI  4 . 7 , 30H  MASS  FLOW ( SLUGS/SQ.FT  .SEC. >  = 
61PE14.7//8H  M  MAX  .  =  1  PE 1 4 . 7 , 9H  DELT AP= 1 PE 1 4. 7 / / / / 1 20H  LOGlOtP/ 

7PT )  LOGIO (RHO/RHOT )  H/HT  Z  A/ASTA 

8R ( SNDSPD )  M  AREA/AREA  STAR) 

DO  15  J  =  1 »  800 

READ  ( 17*7)  PRTL.RHORTL.HRT  »Z  .AR.EM.AREAR 
7  FORMAT  (1P7E17.7) 

RVRSTR=SQRT (AREAR ) 

PRT=10.**PRTL 

RHORT=10.**RHORTL 

Q-l.-HRT 

PUNCH  10.EM.RVRSTR.Q.PRT  * RHORT 
10  FORMAT  ( 5F14.8  ) 

IF  (AREAR-1.0001)  20.20.15 
15  CONTINUE 
20  PRT=PROSTR/PSTAG 
RHORT=RHOSTR/RHOSTG 
EM=  1  • 

RVRSTR= 1 . 0 

Q=ASTAR** 2/(50060. *HS TAG) 

PUNCH  10.  EM, RVRSTR.Q.PRT, RHORT 
STOP 
END 
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APPENDIX  C 

COMPUTER  PROGRAM  FOR  TURBULENT 
BOUNDARY  LAYER  CALCULATIONS 


The  Real  Gas  Turbulent  Boundary  Layer  Program  (RGTBLP)  presented 
here  has  been  developed  for  the  calculation  of  turbulent  boundary 
layers  in  axisymmetric  nozzles  using  air  or  nitrogen  under  conditions 
where  the  gas  no  longer  behaves  like  a  perfect  gas.  The  program, 
valid  for  thick  or  thin  boundary  layers,  is  based  on  the  method 
developed  by  Tetervin  (reference  (.16)  )  and  uses  the  auxiliary 
relations  and  the  air  and  nitrogen  properties  discussed  earlier  in 
this  report.  The  air  properties  are  valid  for  pressuresup  to  100 
atmospheres  and  temperatures  up  to  5000°K.  The  nitrogen  properties 
are  valid  up  to  1000  atmospheres  and  3000°K.  Above  these  values,  the 
program  should  be  used  with  some  caution. 

The  computer  program  has  been  set  up  to  allow  the  calculation  of 
multiple  cases  using  the  same  non-dimensionalized  inviscid  core 
contour.  Thus,  the  gas  properties  and  coordinates  for  the  inviscid 
core  are  read  only  once,  during  the  first  case.  This  eliminates  the 
need  for  tape  switching  and  (if  the  subsonic  inlet  option  is  being 
used)  large  quantities  of  data  cards.  The  supersonic  inviscid  core 
data  is  input  from  the  tape  generated  by  the  Axisymmetric  Core 
Program  (ACP) .  The  inviscid  core  data  for  the  subsonic  inlet  and 
throat  regions  is  input  from  punched  cards.  These  cards  are  required 
only  if  the  subsonic  inlet  is  to  be  calculated.  A  utility  program 
for  punching  these  cards  from  the  isentropic  expansion  tape  used  as 
input  to  the  ACP  program  is  given  in  Appendix  B . 

For  very  low  Mach  numbers  at  the  beginning  of  the  subsonic  inlet, 
it  might  be  necessary  to  generate  additional  data  for  the  subsonic 
inlet  calculation  to  ensure  that  the  interpolation  scheme  used  in  the 
program  does  not  lead  to  large  inaccuracies  near  the  start  of  the 
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inlet.  If  required,  additional  data  at  very  low  Mach  numbers  can  be 
obtained  by  recalculating  the  isentropic  expansion  or  by  interpolating 
the  available  data.  For  the  nozzles  discussed  in  this  report, 
additional  data  was  obtained  by  rerunning  Woolley's  program  with 
finer  temperature  steps  near  the  supply  temperature.  For  the  Mach 
number  20  nozzle,  it  was  also  necessary  to  use  linear  interpolation 
to  generate  even  more  data  between  the  supply  values  and  those  at  the 
lowest  Mach  number  calculated  using  Woolley's  program.  The  total 
number  of  inviscid  core  data  cards  was  85,  163,  and  155  (18  by  linear 
interpolation)  for  the  Mach  number  10,  15,  and  20  nozzles,  respectively. 
Except  for  the  18  points  obtained  by  linear  interpolation  for  the  Mach 
20  case,  the  temperatures  used  in  rerunning  Woolley's  program  to 
generate  the  additional  data  for  the  subsonic  inlet  calculation  have 
been  included  in  Table  3.  Because  of  the  large  number  of  cards 
involved,  the  inviscid  core  data  are  not  included  in  the  listing  of  data 
input  given  in  Table  5. 

In  order  to  maintain  control  over  the  interval  of  integration 
used  throughout  the  calculation  (and  to  prevent  errors  caused  by 
backing  up  toward  the  starting  point  during  the  calculation) ,  the 
program-controlled  variable  interval  of  integration  option  available 
in  the  utility  integration  package  FN0L2  was  not  used.  Instead,  FNOL2 
is  used  in  its  constant  interval  of  integration  mode.  However,  the 
interval  can  be  varied  in  a  practical  sense  by  dividing  the  nozzle 
into  sections  of  arbitrary  length  and  integrating  over  each  successive 
section  with  the  desired  (.constant)  interval  for  that  section.  The 
interval  can  be  increased  or  decreased  in  successive  sections,  as 
preferred.  As  indicated  below,  the  interval  DELX  and  the  X  coordinate 
XEND  (distance  from  the  nozzle  throat)  at  which  the  nozzle  section 
ends  are  read  from  a  data  card.  The  boundary  layer  is  computed  using 
the  interval  DELX  until  X  equals  XEND.  Then,  another  card  with 
values  of  DELX  and  XEND  for  the  next  nozzle  section  is  read  and  the 
process  repeated.  Each  successive  value  of  XEND  must  be  greater  than 
the  previous  one,  so  that  the  calculation  proceeds  continuously  toward 
the  nozzle  exit.  The  last  value  of  XEND  must  be  equal  to  or  greater 
than  the  length  of  the  supersonic  portion  of  the  nozzle  (XEXIT)  in 
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order  that  the  program  will  terminate  properly  (and  proceed  to  the 
next  case,  if  any) .  If  the  yalue  of  XEND  given  is  greater  than 
XEXIT,  then  the  calculation  will  stop  at  XEXIT.  Therefore,  to  simplify 
the  choice  of  the  last  XEND,  it  is  convenient  to  use  some  excessively 
large  value  for  XEND  such  as  999.  (Units  of  XEND  are  feet) . 

PROGRAM  INPUT 

The  data  for  this  program  are  input  from  magnetic  tape  and  from 
punched  cards.  The  tape  input  consists  of  the  inviscid  core  contour 
tape  generated  by  the  Axisymmetric  Core  Program  (ACP) .  This  tape  is 
mounted  on  tape  unit  26.  The  card  input  consists  of  the  following 
six  card  groups,  each  consisting  of  one  or  more  cards. 


Group  1  Required  for  each  case.  This  group  consists  of  one  card 
identifying  the  gas  to  be  used.  Either  AIR  or  NITROGEN 
should  be  punched  beginning  in  column  1.  Since  this  program 
is  set  up  to  calculate  multiple  cases  (if  desired),  it  will 
automatically  return  to  this  READ  after  completing  each  case. 
Therefore,  a  blank  card  is  used  after  the  last  data  set  to 
stop  the  program  normally. 

Group  2  Required  for  each  case.  This  group  is  used  to  label  the 
beginning  of  the  output  with  0-50  lines  of  Hollerith 
information,  as  desired.  Card  1  controls  the  reading  of  the 
Hollerith  information  from  0-50  cards  following  card  1.  If 
multiple  cases  are  being  run,  lines  1  through  LIN  (see  below) 
plus  up  to  10  additional  lines  (.arbitrarily  distributed 
among  the  50  allowed)  used  for  the  previous  case  can  be 
replaced  by  new  lines  input  from  cards.  This  permits 
replacement  of  some  of  the  lines  each  run  without  the  need 
of  re-entering  lines  which  stay  the  same. 

Card  Card  Format  Variable  „ .  .  . 

Number  Columns  Used  Name  e  inl  10n 


1 


1-5  15  LIN  Number  of  cards  (lines) 

of  labeling  to  be  read  in. 
Do  not  include  replace¬ 
ment  cards  covered  by 
LREPLC  (see  below) . 

6-10  15  LOUT  Total  number  of  lines  of 

labeling  to  be  written 
out . 
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Card 

Number 

Card 

Columns 

Format 

Used 

Variable 

Name 

Definition 

11-60 

1015 

LREPLC 

Line  numbers  of  0-10 
lines  from  previous  case 
to  be  replaced  for  this 
case. 

2-51 

1-72 

12A6 

TITLE 

Hollerith  labeling,  72 
characters  per  card. 

Number  of  cards  (0—50) 
equal  to  LIN  plus  one 
for  each  non-zero  value 
read  for  LREPLC. 

Required  for  each  case.  Two  cards  in  this  group  contain  most 
of  the  required  initial  values  and  control  numbers.  Several 
of  the  control  numbers  have  default  values,  such  that  if  the 
value  read  from  the  card  is  zero,  a  particular  value  is 
assigned  in  the  program..  This  eliminates  the  need  to  punch 
control  numbers  for  the  most  used  options. 

Card 

Number 

Card 

Columns 

Format 

Used 

Variable 

Name 

Definition 

1 

1-10 

F10. 5 

THETA 

Momentum  thickness  (.ft) 
at  initial  station 

11-20 

F10. 5 

REXIT 

Radius  (ft)  of  inviscid 
core  at  nozzle  exit 

21-30 

F10 . 5 

XPUNCH 

Axial  distance  (ft)  from 
throat  at  which  punching 
of  nozzle  coordinates  intc 
cards  is  to  begin 

31-40 

F10. 5 

XO 

Axial  distance  (ft)  from 
throat  at  which  calcula¬ 
tion  begins.  XO<0  for 
subsonic  inlet  calculation 

41-50 

F10. 5 

RO 

Radius  (ft)  of  subsonic 
inlet  at  XO.  Ignored  if 
XOhO. 

51-60 

F10. 5 

OMEGAO 

Angle  (degrees)  of  sub¬ 
sonic  inlet  at  X0.  OMEGAO 
<0  for  converging  wall  at 
XO.  Ignored  if  XO>0. 

61-70 

F10 . 5 

III 

Shape  parameter  at  initial 

station.  If  HI=0,  program 
calculates  starting  value 
using  THETA. 
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Card  Card  Format  Variable 

Number  Columns  Used  Name 


Definition 


1-5 

15 

NN 

Identification  number. 

6-10 

15 

LM 

Tape  output  control.  If 
LM>0,  output  will  be  on 
tape  16  as  well  as  on 
system  tape. 

11-15 

15 

NFS 

Number  of  files  to  be 
skipped  before  reading 
inviscid  core  contour 
data  from  tape  26. 

16-20 

15 

IP 

Overridden  by  Group  6. 

Can  be  left  blank. 

21-25 

15 

J 

Integration  mode  control. 
If  1,  Runga-Kutta  only. 

If  2,  Runga-Kutta  for 
first  four  and  last 
intervals,  Adams-Moulton 
for  others.  If  3,  Runga- 
Kutta  repeated  for  each 
step.  (See  reference  (26) 
for  details) .  The  Adams- 
Moulton  mode  ( J=2)  is 
normally  used. 

26-30 

15 

LAWCF 

Skin  friction  law  control 
If  1,  Spalding-Chi  lav/. 

If  2,  Ludwieg-Tillrrtann 
lav/  with  reference 
enthalpy . 

Group  4  Required  only  for  the  first  case  and  even  then  only  if  the 
subsonic  inlet  option  is  to  be  used  in  one  of  the  cases. 
This  group  of  cards  contains  the  inviscid  core  properties 
required  for  the  subsonic  inlet  calculation.  The  form  used 
for  these  properties  is  the  same  non-dimensionalized  form 
used  in  the  Ax i symmetric  Core  Program  (ACP)  for  the  super¬ 
sonic  portion  of  the  inviscid  core.  A  utility  routine  for 
obtaining  these  cards  from  the  .isentropic  expansion  tape 
used  as  input  to  the  ACP  program  is  given  in  Appendix  B. 

The  last  card  of  this  group  must  correspond  to  the  throat 
conditions,  i.e.  the  Mach  number  must  be  one. 
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Group  5 


Card  Card  Format  Variable 
Number  Columns  Used  Name 


Definition 


EACH  1-14 


15-28 


29-42 


43-56 

57-70 


FI 4, 8  EMW 

F14.8  RyRSTR 


F14.8 


F14.8 


QW 


PW 


FI 4. 8  RORSTG 


Mach,  number 

Radius/ throat  radius 
Ci.e.  square  root  of 
area  ratio) 

Ratio  of  square  of  flow 
velocity  to  twice  the 
supply  enthalpy 

Ratio  of  pressure  to 
supply  pressure 

Ratio  of  density  to 
supply  density 


Required  for  each  case.  This  card  defines  the  wall  tempera¬ 
ture  distribution  to  be  used  in  the  calculation.  The  three 
options  available  are:  constant  temperature  wall,  adiabatic 
wall,  and  prescribed  temperature  distribution  based  on 
equation  (25).  Variables  TX(3),  TX(4),  TX(5),  TX(6),  TX(8) 
are  ignored  if  prescribed  temperature  distribution  is  not 
used . 


Card  Card  Format  Variable 

Number  Columns  Used  Name 


Definition 


1 


1-5 


6-15 


16-25 


26-35 


15  M  Wall  temperature  option 

control.  If  1,  constant 
wall  temperature.  If  2, 
prescribed  wall  tempera¬ 
ture  distribution.  If  3, 
adiabatic  wall  conditions 


F10.5  TX(2)  Constant  wall  temperature 

( °K)  if  M= 1 .  Throat  wall 
temperature  (°K)  if  M=2. 
Starting  guess  for  wall 
temperature  (°K)  at  initial 
station  if  M=3. 


F10.5  TX(3)  Distance  (ft)  from  throat 

of  first  intermediate 
temperature  point  for 
prescribed  distribution. 

F10.5  TX(4)  Wall  temperature  C°K)  at 

TX  (3 )  . 
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Card  Card  Format  Variable 

Number  Columns  Used  Name 


Definition 


36-45  P10.5  TX(5) 


46-55  F10.5  TX(6) 


Distance  (ft)  from  throat 
of  second  intermediate 
temperature  point  for 
prescribed  distribution. 

Wall  temperature  (°K)  at 
TX  ( 5 )  . 


56-65  F10.5  TX ( 8 ) 


Wall  temperature  (°K)  at 
nozzle  exit. 


Group  6  Required  for  each  case.  These  cards  control  the  interval  of 
integration  (see  discussion  above)  and  the  frequency  of 
output  for  each  section  of  the  nozzle. 


Card  Card  Format  Variable 

Number  Columns  Used  Name 


Definition 


EACH  1-10  FlO. 5  DELX  Interval  of  integration 

(ft). 

11-20  FlO. 5  XEND  Distance  (ft)  from  throat 

to  downstream  end  of  section 
in  which  DELX  will  be 
used.  If  XEND  is  greater 
than  length  of  supersonic 
portion  of  nozzle  (XEXIT) , 
calculation  will  terminate 
at  XEXIT. 

21-25  15  IP  Output  v/rite  frequency, 

e.g.  output  at  every  other 
step  if  IP=2.  Regardless 
of  value  of  IP,  output  will 
be  written  for  first  and 
last  steps.  If  zero, 
program  sets  IP=1. 


PROGRAM  OUTPUT 


Output  written 
identification  from 


on  the  system,  tape  will  include  all  pertinent 
the  ax i symrae.tr ic  core  contour  tape  used  as  input 
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plus  nearly  all  the  card  input  quantities  listed  in  the  previous 
section.  In  addition,  for  each  step  of  the  integration  at  which 
output  is  written  Csee  IP  in  input  card  Group  6) ,  a  table  of  forty- 
eight  quantities  is  written.  This  table  includes  the  gas  properties 
at  the  edge  of  the  boundary  layer,  the  various  boundary  layer  thick¬ 
nesses,  heat  transfer  parameters  and  other  quantities  of  interest. 

Each  page  of  tables  is  labeled  with  headings  for  easy  identification 
of  the  various  quantities  and  the  headings  are  defined  at  the 
beginning  of  the  output.*  If  a  permanent  copy  of  the  output  is 
desired,  the  input  quantity  LM  should  be  made  greater  than  zero.  In 
this  case,  output  identical  to  that  written  on  the  system  tape  will 
be  written  on  tape  16.  In  addition  to  the  output  written  on  tape, 
the  user  can  obtain  the  wall  and  core  coordinates  on  punched  cards. 

This  card  output,  governed  by  the  two  card  input  quantities  IP  and 
XPUNCH,  consists  of  the  five  quantities  X,  RW,  R,  XINCK,  RWINCII,  defined 
as  follows: 

X 

RW 

R 

XINCH 

RWINCH 


distance  (ft)  from  nozzle  throat  along  centerline 
radius  (ft)  of  nozzle  wall  at  X 
radius  (ft)  of  inviscid  core  at  X 
X  expressed  in  inches 
RW  expressed  in  inches 


*The  reader  can  identify  the  forty-eight  quantities  tabulated  by 
comparing  the-  WRITE  statement  in  subroutine  OUT  with  the  list  of 
definitions  contained  in  subroutine  LABEL3. 
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PROGRAM  LISTING 


C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT»XT,HO, 

1  DXDX.DRDX ,DEMDX,DRHODX»DPEDX,DUEDX,DHEDX, 

2  CF»COEFF»DELSTR, DELTA, DHIDX,DH1 » DH2  » DRWDX  » 

3  DSVDEL,DSVTH,DSWDX,DTHDX,EMUE,GAS,HAD,HFI * 

4  HI » HJ  »HK  *HW  , IP,IP1,LAWCF,LM»M,N»PAGE,PR»Q, 

5  re,recvry,*reth,reyanl>rhow,roe,rw,st  * SUM1  * 

6  SUM2,TE,THETA,THVDEL,TW,TX(8) ,  XEND , XP UNCH , 

7  XTRA1*XTRA2.XTRA3*ZE 
DIMENSION  Z(5) »D( 5) 

DATA  A  I R/6HA I R  /»  GAS2  / 6HN I TROG/ » BLANK / 6H 

C 

C  PRANDTL  NUMBER  ASSUMED  CONSTANT  AT  .72 

PR  =  «  72 

C  RECVRY  =  PRANDTL  NUMBER  TO  1/3  POWER 

RECVRY=. 89628 

C  REYANL  =  PRANDTL  NUMBER  TO  -2/3  POWER 

REY  ANL= 1 • 2448 
C 

10  READ  (5,20)  GAS 
20  FORMAT  (A6) 

IF  (GAS. EQ. BLANK)  GO  TO  100 
IF  (GAS. EQ. AIR)  GO  TO  40 
IF  ( GAS.EQ.GAS2  )  GO  TO  40 
WRITE  (6,30)  GAS 

30  FORMAT  ( 1H1  A6.22H  IS  NOT  AN  ALLOWED  GAS  ) 

GO  TO  100 
C 

40  CALL  LABEL1  (GAS) 

C 

READ  ( 5 ,50  > THET A, REX  IT, XPUNCH.XO.RO, OMEGA 0»H I ,NN  »  LM  .NFS  » I 
50  FORMAT  (7F10. 5/615) 

IF  (J.EQ.O)  J=2 
IF  (IP.EQ.O)  IP=1 
IF  (LAWCF.EQ.O)  LAWCF=1 
IF  (NFS.EQ.O)  GO  TO  70 
DO  60  N F  =  1 ,NFS 
EXTERNAL  UN26 
CALL  FSFILE  (UN26) 

60  CONTINUE 
70  X=XO 
R  =  R0 

DRDX  =  T  AN  (OMEGAO/57. 29578 ) 

I P 1 = 1- I P 
PAGE=0 • 

C 

CALL  CORE  ( X , LM  ,  1  ) 

M  =  4 

CALL  WALL  (X) 

CALL  LABEL2  ( X 0 , RO , OMEGAO , NN , J ) 

CALL  CORE  ( X,LM,2  ) 

CALL  LABEL3  (LM) 

C 

IF  (HI.GT.O.  )  GO  TO  80 
CALL  RHOTZE 

IF  (GAS. EQ. AIR)  CALL  EMUAIR  ! HE , PE , EMUE ) 

IF  (GAS. NE. AIR)  CALL  EMUN2  ( T E  ,  RHOE  ,  EMUE ) 

RETH=UE*RHOE*THETA/EMUE 

RETHLG=ALOG10 ( RETH ) 

HFI=1Q.**(0.599-RETHLG*(0«198-0.0189*RETHLG) ) 
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H I =HF  f 

80  Z ( 1 ) =THET A 
Z ( 2 ) =H I 
CF  =  0. 

90  READ  (5,91)  DELX,XEND,IP 

91  FORMAT  <2F10.5,I5) 

IF  (IP.EQ.O)  IP=1 
EXTERNAL  DER, TERM, OUT 

CALL  FN0L2  ( J , 2 , DELX  ,  0  ,  I P  ,  0 , X ,Z , D , DER , TERM , OUT ) 
IF  ( XEND.LT.XEXIT )  GO  TO  90 

95  END  FILE  16 
GO  TO  10 

100  REWIND  16 
REWIND  26 
STOP 
END 


FUNCTION  COMP ( U , X ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  EVALUATES  Z  FROM  LOGlO  P  AND  LOGlO  RHO 
AND  IS  IDENTICAL  TO  COMP  USED  IN  AIETP  DATED  7/01/65 

I F ( U-0 .5)10,5,5 
5  Y=-l • 0 
GO  TO  15 

10  Y=1,0 

15  A0=1.042+U*( 1.069  +  U*! 6.3 54E-3+U*( —3.3 E-3+U*( —9.15 5E-4+U*<  5.952 E-4 
l+U*(2.866E-4+U*(4.287E-5+2.18E-6*U)  )))))) 

Al=-44.71+U*( 8 . 358+U* ( - . 2 31 8+U* ( -. 1 07  +  U* ( - . 1 1 76+U* ( 1.642E-2+U*< 

1 2  •  46 1  E  - 2  +  U#  (  5  •  2 1  E  - 3  +  3  •  362  E  — 4*U  ) ')  )  )  )  )  ) 

CO=l . 137+U*< 1.063+U*( 2.836 E-2+U*( 1 . 698E-2+U* ( -2 . 2 3E-4+U* ( -2 . 836E-3 
l+U*(-8.384E-4+U*< -9 . 2 59 E-5-3 . 45 1 E-6*U )  )))))) 

C2=.1805+U* (-9.217E-2+U*( 1 . 07 5E-2+U* ( 3 ,633 E-3+U* ( -1 . 602E-4+U* ( 
l-4.084E-4+U*(-9.916E-5  +  U*(-8.743E-6-2.119E-7*U)  )))))) 

EO  =  l ,43  5  +  U* ( 1. 1  +  U*( 3.016E-2+U*( 1 . 447E-3+U* ( -6 . 107 E- 3+U* ( -2 . 075 E-3 
l+U*(-1.37E-4+U*(2.689E-5+3.088E-6#U> )))))) 

ZO=-l  5  •  *■  (  U-  •  5  ) 

IF(ZO-10. >21,20,20 

20  El =-22 «  0 
GO  TO  30 

21  IF(ZO+8. 122,22,23 

22  El=-56,6+p,8*U 
GO  TO  30 

23  El=-22.0+(9.8*U-34.6) / ( 1 ,+EXP  ( ZO ) > 

30  E2=-.2476+U*(-4.8  58E-2+U*( 1.217E~2+U*<  3 . 934E-3+U* ( -3 . 70  5E-4+U* ( 
l-5.415E-4+U*(-1.095E-4+U*(-5.681E-6+1.694E-7*U) )))))) 

E3= . 3904+U* ( .3021+U*(4.053E-2+U*(-1.283E-2+U*(-1.935E-3+U*( 
15.477E-4+U*(2.667E-4+U*(4.022E-5+2.142E-6*U) )))))) 

G0=1.87016+U*< 1 ,08549+U*( 2 . 39061E-2+U* ( 2 . 06972E-2+U* ( 1.71414E-2+U 
l*(6.45427E-3+U*( 1 . 0 58 3 5 E-3+6 . 2 73 06 E-5*U ) ) ) ) ) ) 

Gl=-13.9930+U*( 2.33941+U* ( .431717+U*( . 162394+U* ( 5 .00809E-2+U* ( 

17, 94148E-3  +  4. 80406E->4*U  )  )  )  )  ) 

G2  =  . 408  95  5  +  U* ( -1 . 8 3 20 5 E-2  +  U* ( -9 . 9 2 1 56E-3  +  U* ( 5 . 947 58 E-3  +  U* ( 
14.30655E-3+U*(9.39536E-4+U*(8.33593E-5+2.33120E-6*U) ) ) ) ) ) 
G3=-.718311+U* (-.423658+U* ( 2 . 64678 E-2+U* ( 1 . 1 7670 E-2+U* ( 
11.67085E-3+U*( 5 . 8 5089E-5-2 . 2 1 499E— 6*U ) ) ) ) ) 

F I  0  =  2  «  462  +  U* ( 1.03 5+U* ( -7.0  8  3E-3  +  U* ( 6 . 8 93E-3+U* ( 3 . 469E-3  +  U* ( 
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l-1.743E-4+U*(-4.025E-4+U*(-8.534E-5-5.571E-6*U) )))))) 

FI 1=-16.37+U* ( 1. 56 1+U* ( -• 2057+U* ( -5  «  55  8E-2+U* ( 4. 067E-2+U* ( 1 . 99  3E-2 
l+U*(-2.783E-3+U*(-1.926E-3-1.952E-4*U)  )))))) 

F I  2  = .  1 9 86  +  U* ( -8 . 448 E-2+U* ( -6 . 44E-4+U# ( 2 . 533E-4+U* ( -6. 545E-4+U* ( 
l-5.846E-5  +  U*(8.10  9E-5  +  U*<  2.08 5E-5  + 1.431 E-6*U) )))))) 

40  Z 1 = A 1* ( X-AO ) 

IF(Z1-10. 142,41,41 

41  F1  =  0.0 
GO  TO  51 

42  IFIZ1+8. >43*43*44 

43  Fl=.3468*( X-AO) 

GO  TO  50 

44  I F ( ABS  (Zl)-. 01)45*45*46 

45  Fl=-.3468/Al 
GO  TO  51 

46  Fl=.3468*(X-AO>/< l.-EXP  (Zl)) 

50  Z2=-100.*< X-CO) 

I F ( Z2-10. 152*51*51 

51  F2  =  0 . 0 
GO  TO  61 

52  IF(Z2+8. 153,53*54 

53  F2=C2* ( X-CO ) 

GO  TO  60 

54  I F ( ABS  (Z?)-. 01 155.55,56 

55  F2=.01*C2 
GO  TO  61 

56  F2=C2*(X-CO)/(l.-EXP  (Z2)) 

60  Z3=E1* ( X-EO ) 

IF ( Z3-10. ) 62*61 ,61 

61  F3=0.0 
GO  TO  81 

62  IF(Z3+8. 163,63,64 

63  F3=E2*X+E3 
GO  TO  80 

64  I F ( Y )  65,70,70 

65  IFtABS  <Z3)— .01)  66,66,67 

66  F3=-E2/E1 
GO  TO  81 

67  F3  =  E2*(X-EO)/( l.-EXP  <  Z3 )  > 

GO  TO  80 

70  F3= ( E2*X+E3 ) / ( l.+EXP  (Z3)> 

80  Z4=G 1* ( X-GO ) 

I  F ( Z4- 10. >82,81,81 

81  F4  =  0 . 0 
GO  TO  91 

82  IFIZ4+8. 183,83,84 

83  F4=G2*X+G3 
GO  TO  90 

84  F4= ( G2#X+G3 ) / ( 1 .  +  EXP  ( Z4  )  ) 

90  Z  5  =  F 1 1* ( X-F I  0 ) 

IFIZ5-10. 192,91,91 

91  F  5  =  0 . 0 

GO  TO  100 

92  IFIZ5+8. 193 ,93,94 

93  F5  =  F I  2# ( X-F I  0 1 
GO  TO  100 

94  I F ( ABS  ( Z  5 1  - . 01 195,95,96 

95  F5  =  -FI 2/FI  1 
GO  TO  100 

96  F5  =  FI 2* ( X-FIO ) / ( 1 .-EXP  ( Z 5 1  1 
100  RECIP=1,-F1-F2+F3+F4+F5 

COMP= 1 . /REC I P 
RETURN 
END 
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SUBROUTINE  CORE  ( X , LM  *  L 1  ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.6L0WACKI 

COMMON  QP*REXIT  » X  EX  I  T  »XTRE»HSTAG*DQPDX 

DIMENSION  Q  (  7  )  »  01(7),  Q2 ( 7 ) ,  Q3 ( 7 )  ,  OP ( 7 )  ,  DQPDX(7),  A ( 7  > ,  B(7), 
1KP ( 999 )  »XW(999) ,YW(999)  ,EMW(999)  ,  PW ( 999 )  ,QW ( 999 )  ,R0RSTG(999) 

DATA  NCASE/0/ 


Q (  1  )  =X  ,  Q ( 2  )  =R  ,  Q ( 3 )  =  EM  ,  Q(4)=RH0E,  Q(5)=PE,  Q(6)=UE,  Q(7)=HE 


GO  TO  ( 10,200,390)  ,  LI 

READ  CORE  PROPERTIES  DURING  FIRST  CASE  ONLY 
10  NCASE=NCASE+1 

IF  (NCASE.GT.l )  GO  TO  40 

READ  HEADINGS  IDENTIFYING  EXPANSION  AND  CORE  PROGRAM  CALCULATIONS 
READ! 26 , 2 0 ) NN  ,  PST AG , T T , RHOSTG , HST AG ♦ T Z , SST AG , T AR , TST AR , AST AR , 
1PR0STR,RH0STR,RH0QST,EMMAX,DELTAP 

20  FORMAT ( 37H1  ISENTROPIC  EXPANSION  TABLE  NUMBER  I4//18H  SUPPLY  COND 
1ITI0NS/9H  P(ATM. ) = 1  PEI 4 . 7 , 1 1 H  T ( KEL V I N )  =  1 P E 1 4 . 7 , 2 2H  RHO ( SLUGS/CUB  I 
2C  FT. )=1PE14.7/23H  ENTHALPY ( B TU/ LB . MASS ) = 1 PE 1 4 . 7 , 4H  Z=1PE14.7, 

36H  S/R  =  1PE14.7,10H  A / AS T AR= 1  PE  1 4 . 7/ / 1 8H  THROAT  COND I T I ONS / 1 1H  T( 
4KELVIN ) =1PE14.7  »18H  ASTAR(FT./SEC.)=1PE14.7,10H  P(ATM.)=1PE14.7/ 

522H  RHO (SLUGS/CUBIC  FT . ) = IP E 1 4 . 7 , 30H  MASS  FLOW (SLUGS/SQ.FT .SEC. )- 
61PE14.7//8H  M  MAX . = 1 P E 1 4. 7 , 9H  DELT AP= 1 P E 1 4 . 7/ / / / 1 20H  ****** 

7********************************* 

8*  **************  7*  ****** 

READf  26,30 ) NNN , EMT , RSTAR , XT, RCT, EMZP, I FTH , P A , PB , P C , PD , KT  ,  KK , GAM A » 
1XZ,DZX,DTX,EPSA,EPSB 

30  FORMAT ( 46H0AXIALLY  SYMMETRIC  NOZZLE  WALL  CONTOUR  NUMBER  I3///23H  T 
1ERMINAL  MACH  NUMBER  =F10.6»16H  THROAT  RADIUS  =1PE14.7//5H  XT  =1PE1 
24.7  »10H  RCT  =1PE14.7,11H  EMZP  =1PE14.7//6H  IFTH=I2*7H  P 

3 A  =  1 PE14 . 7 , 7H  PB= 1  PE  1 4. 7 , 7H  PC  =  1 PE 1 4 . 7 , 7H  PD  =  1PE14.7 ,7H 

4  KT=I3,4X,3HKK=I3//6H  GAM A= IP E 1 0 . 3 , 4X , 3HXZ= 1 P E 1 0 . 3 , 4X , 4HDZX = IP E 1 0 . 
53»4X»4HDTX=1PE10.3»4X,5HEPSA=1PE10.3*4X»5HEPSB=1PE10«3/13H1  NN  K 
6  XW  1 4X  » 2  H  YW  14X,6HTHETAW  10X,2HMW  14X,9HPW/-P  STAG  7X.2HQW  14X, 
712HRHO/RHO  STAG) 

C 

40  WRITE(6,20)NN,PSTAG,TT,RHOSTG,HSTAG,TZ»SSTAG»TAR*TSTAR»ASTAR* 
1PROSTR , RHOSTR , RHOQST , EMMA X , DELT AP 
WRITE ( 6,50 )NNN,EMT, RSTAR, XT, RCT, EMZP, I FTH , P A , PB , P C , PD , KT , KK  , GAMA , 
1XZ»DZX,DTX,EPSA ,EPSB 

50  FORMATf 46H0AXI ALLY  SYMMETRIC  NOZZLE  WALL  CONTOUR  NUMBER  I3///23H  T 
1ERMINAL  MACH  NUMBER  =F10.6,16H  THROAT  RADIUS  =1PE14.7//5H  XT  =1PE1 
24.7, 10H  RCT  =1PE14.7,11H  EMZP  =1PE14.7//6H  IFTH=I2,7H  P 

3A= 1 PE14 . 7 , 7H  P B=  1 P E 1 4 . 7 , 7H  P C= 1 P FI  4 . 7 , 7H  PD= 1 P E 1 4 . 7 , 7H 

4  KT=I3,4X,3HKK=I3//6H  GAMA  =  IP  El  0 . 3 » 4X , 3HXZ  =  1  PE  10 . 3  * 4X , 4HDZX  =  1  PE  10 . 
53,4X,4HDTX=1PE10.3,4X,5HEPSA=1PE10.3,4X,5HEPSB=1PE10.3/120H  *  *  * 
6********************************* 
7************************) 

IF  (LM.EQ.O)  GO  TO  60 

WR I TE ( 16,20)NN,PSTAG,TT,RH0STG,HSTAG,TZ,SSTAG,TAR,TSTAR,ASTAR, 
1PROSTR , RHOSTR , RHOQST, EMM AX,  DEL TAP 
WRITE! 16,50 ) NNN, EMT, RSTAR,  XT, RCT, EMZP  ,IFTH,PA,PB,PC,PD,KT,KK, 

1GAMA  ,XZ  ,DZX  ,DTX,EPS’A  ,EPSB 
60  IF  (NCASE.GT.l)  GO  TO  120 
C 

STAGH=50060.*HSTAG 
STAGP=2116.*PSTAG 
XW ( 1 )=0. 

YW ( 1 ) =R  STAR 
EMW ( 1 >=1. 

PW( 1 )=PROSTR/PSTAG 
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QW ( 1 )=ASTAR**2/STAGH 
RORSTG ( 1 ) =RHOSTR/RHOSTG 
1=1 
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C 

IF  (X.GE.O. ) GO  TO  100 

C  READ  INLET  CORE  PROPERTIES  FROM  CARDS 

1=0 

70  1=1+1 

READ  (5,80)  EMW(  I  ) ,RVRSTR,QW(  I  ) *PW(  I  )  >RORSTG( I ) 

80  FORMAT  (5F14.8) 

YW ( I ) =RVRSTR#RSTAR 
IF  ( EMW (  I ) .LT.l.  )  GO  TO  70 
12=1 

READ  CHARACTERISTIC  CORE  CONTOUR  AND  PROPERTIES  FROM  TAPE 
100  IF  (I.EQ.999)  GO  TO  500 
1  =  1  +  1 

READ(26,110)NN,KP( I  ) , XW ( I )  ,  YW (  I  )  , THET  AW  ,  EMW (  I  )  ,  PW ( I )  ,  QW (  I  )  , 

1RORSTG ( I ) 

110  FORMAT  ( 2  I  4  »  7E 16  «  7 ) 

IF  (KP( I ) .LT.KT)  GO  TO  100 
XWL A  ST  =  XW ( I  ) 

120  XEXIT=XWLAST*REXIT 
XTRE  =  XT*REX  I  T 
J  =  0 
L  =  1 

IF  (X.GE.O.  )  GO  TO  300 

INTERP  WILL  INTERPOLATE  TO  FIND  VALUES  AT  BEGINNING  OF  INLET 
130  XWO=X/REX I T 

YW0=QP(2)/REXIT 
DYDX0=DQPDX(2) 

JMAX=999 

CALL  INTERP  ( 3  ,  JMAX , YW , EMW , YWO , EMWO » J ) 

CALL  INTERP  < 3  ,  JMAX , YW ,QW , YWO » GWO » J 1 ) 

CALL  INTERP  ! 3  ,  JMAX , YW , PW , YWO * PWO * J 1 ) 

CALL  INTERP  < 3 , JMAX > YW , RORSTG » YWO * RORS TO » J 1  ) 

SUBROUTINE  SHAPE  WILL  DETERMINE  SHAPE  OF  INLET  CORE  CONTOUR 
CALL  SHAPE  ( XWO , Y WO , DYDXO , XW , YW , JMAX , J , I  2 ) 

START  INTERPOLATIONS  AT  BEGINNING  OF  INLET 
Q ( 1  )  =REXIT*XWO 
Q(  2)=REXIT*YW0 
0(3) =EMW0 

Q( 4) =RHOSTG*RORSTO 
Q(5)=STAGP*PW0 
Q ( 6 ) =  SQRT ( STAGH*QW0  ) 

Q ( 7 ) =HST  AG* ( l.O-OWO ) 

J  =  J  +  1 
GO  TO  310 

WRITE  INLET  CORE  CONTOUR  AND  PROPERTIES  INPUT  FROM  CARDS 
200  IF (X.GE.O.)  GO  TO  600 
WRITE(6,210) 

210  FORMAT  (1H1  13X.85H SUBSONIC  INVISCID  CORE  DATA  REQUIRED  FOR  CALCUL 
* AT  I NG  BOUNDARY  LAYER  IN  SUBSONIC  INLET  /1H0  13X,85H(FORM  USED  IS  S 
*AME  NON-DIMENSIONALIZED  FORM  USED  IN  ACP  FOR  SUPERSONIC  INVISCID  C 
*ORE)  / / 1 HO  8X , 1H I  8X.2HXW  14X.2HYW  14X.2HMW  11X.9HPW/P  STAG  10X, 

1  2HQW  9X.13HRHOW/RHO  STAG) 

IF  (LM.GT.O)  WRITE  (16,210) 

DO  230  11=1,12 

WRITE  (6,220)  I  1  ,  XW (  II ) »YW(  II)  ,  EMW (  I  1 )  , PW (  1 1  )  ,  QW !  I  1  )  ,  RORSTG (Ill 
220  FORMAT  (  1 1 0 , 1 P6E 1 6 . 7 ) 

IF  (LM.EQ.O)  GO  TO  230 

WRITE  (  16,2  20 )  I  1 ,XW (  I  1 ) »YW(  1 1 )  ,EMW<  I  1 )  , P W (  I  1 ) , QW (  I  1 ) , RORSTG (  1 1 ) 
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230  CONTINUE 
GO  TO  600 
C 

300  J=J+1 

Q< 1 )=REXIT*XW( J) 

Q( 2 )=R£XIT*YW(J) 

Q( 3 ) =EMW ( J ) 

Q(4)=RH0STG*R0RSTG< J) 

Q(5)=STAGP*PW(J) 

Q ( 6 ) =SQR  T  ( STAGH*QW ( J )  ) 

Q ( 7  )  =HSTAG* ( 1 . 0-QW ( J )  ) 

GO  TO  (310,330,350) ,L 
310  DO  320  1=1,7 
320  Q1 ( I ) =Q ( I  ) 

1  =  2 

GO  TO  300 
330  DO  340  1=1,7 
340  Q2( I )=Q( I ) 

L  =  3 

GO  TO  300 
350  DO  360  1=1,7 
360  Q3 ( I ) =Q (  I  ) 

C  QP( I >=A< I )*R**2+B( I )*R+Q1( I ) 

370  A 1 =Q2 ( 1 )-Q3 ( 1  ) 

A2  =  Q1 ( 1 ) -Q3 ( 1  ) 

A3  =  Q1( 1 ) -Q2 (  1  ) 

DEN  =  A 1* A2*A  3 
DO  380  1=2,7 

A< I )  =  ( ( Q3 (  I >  — Q 1 ( I ) ) *A3-(Q2 ( I ) -Q1 ( I ) ) *A2 ) /DEN 
380  B ( I ) = ( Q 1 ( I > -Q2 ( I ) )/A3+A( I )*A3 
390  QP( 1 ) =X 

I F ( QP ( 1 ) -Q3 ( 1  )  )  400,400,420 
400  R  =  QP ( 1 ) -Q1 ( 1  ) 

DO  410  1=2,7 

QP (I) =Q1 ( I ) +  R* ( B (  I ) +  R*A ( I  )  ) 

410  DQPDX( I )=2.0*A< I ) *R+B ( I  ) 

GO  TO  600 

420  IF(KT.LE.KP( J) )  GO  TO  400 
J  =  J+1 

Q ( 1 ) =R  EX  I T*XW ( J ) 

Q(2)=REXIT*YW( J) 

Q ( 3 ) =EMW ( J ) 

Q( 4) =RHOSTG*RORSTG( J) 

Q(5)=STAGP*PW(J) 

Q ( 6 ) =SQR  T  ( STAGH*QW ( J )  ) 

Q ( 7 ) =HS  T  AG* ( 1 • 0-QW ( J)  ) 

DO  430  1=1,7 
01 (  I ) =02 ( I  ) 

Q2 (  I  ) =Q3 (  I  ) 

430  Q3( I )=Q( I ) 

IFtQP (11-03(1) 1370,370,420 
500  WR ITE(6»510) 

510  FORMAT  ( 54H  STORAGE  REQUIRED  BY  INPUT  DATA  EXCEEDS  THAT  ASSIGNED) 
600  RETURN 
END 


SUBROUTINF  DER(X,Z*D> 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT,XT,HO, 

1  DXDX,DRDX,DEMDX,DRHODX,DPEDX,D'JEDX,DHEDX, 
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2  CF,C0EFF,DELSTR»DELTA,DHIDX,DH1 ,DH2 ,DRWDX> 

3  DSVDFL»DSVTH*DSWDX»DTHDX»EMUE»GAS»HAD»HFI  . 

4  HI»HJ*HK»HW»IP»IP1»LAWCF»LM»M*N»PAGE»PR»Q» 

5  RE»RECVRY»RETH.REYANL,RH0W,R0E,RW,ST,SUM1 » 

6  SUM2»TE.THETA.THVDEL,TW,TX(8) » XEND , XP UNCH , 

7  XTRA1 »XTRA2»XTRA3»ZE 

8  ,CDLSTR,DLSTRR*THETAR,DTHRDX 
DIMENSION  2(5)  , D ( 5  > 

DATA  XPREV/-100 . / »  AIR/6HAIR  /,  APREV/-1 . / 

C 

C 

C  COMPUTE  FUNCTIONS  OF  X  ON  FIRST  PASS  ONLY 

10  IF  (X.EQ.XPREV)  GO  TO  20 
XPR  E V  =  X 
Z ( 3 ) =X 

CALL  CORE  ( X  *  0  *  3 ) 

RHOEUE=RHOE*UE 
DLNUE=DUEDX/UE 
DLNRHO=DRHODX/RHOE 
HEVH0=HE/H0 
DH2=RECVRY* ( HO-HE ) 

HAD=HE+DH2 
CALL  WALL  (X) 

DH1=HAD-HW 

FOFX=( 0.5*DHEDX/HE-DLNUE) *HEVH0 

GOFX=32.174*RHOFUE*DHl 

CALL  RHOTZE 

IF  (GAS. EQ. AIR)  CALL  EMUAIR  ( HE  >  PE  »  EMUE ) 

IF  (GAS. NE. AIR)  CALL  EMUN2  ( T E  ,  RHOE , EMUE ) 
RE=RHOEUE/EMUE 
N=  1 
C 

C  COMPUTE  FUNCTIONS  DEPENDING  ON  THETA  AND  HI 

20  THETAR=Z( 1 ) 

H  I  =  Z  (  2  ) 

HJ=HI+1 . 

HK= (HI-1.1/2. 

CALL  GAUSS 
A  =  0  • 

DO  30  1=1 » 20 

IF  ( THETA. NE.THETAR  )  A  =  DE LT A/ R W/ DSWDX 

THRVDL=THVDEL-A*( 2 . *X TR A 2 ZH J-X TR A3 / H I ) 

DELTA=THETAR/THRVDL 

THETA=THVDEL*DELTA 

DELSTR=DSVDEL*DELTA 

DSRVDL=DSVDEL-A*( . 5-2 . *XT R A2 / H J ) 

DLSTRR=DSRVDL*DELTA 

CDLSTR=RW*DSWDX*( 1. -SORT ( 1.-2. *A*DSRVDL) ) 

IF  (CDLSTR.LE.O. )  CDL ST R= DELS T R 
CALL  SLOPE  (X.CDLSTR.DRDX.DRWDX) 

DSWDX=SQRT ( 1 ,+DRWDX**2 ) 

RW=R+CDLSTR /DSWDX 
CONTRL= 1 . 

IF  (A  .NE.  0.)  CONTRL=ABS< 1 ,-APREV/A) 

IF  (CONTRL  .LE. .0001. AND. I.GT.2)  GO  TO  40 
APREV=A 
30  CONTINUE 
40  RETH  =  RE*THET A 

RETHLG=ALOG10( RETH) 

HFI=10.**(0.599-RETHLG*(0.198-0.0189*RETHLG) ) 

CALL  SKIN 
HALFCF=CF/2. 

ST=HALFCF*REYANL 
Q=ST*GOFX 
N  =  2 
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C  EVALUATE  required  derivatives  of  thetar  and  hi 

DTHRDX=HALFCF*DSWDX-THETAR*(DLNUE*(2.+DSVTH)+DLNRH0+DRWDX/RW) 
DHIDX=F0FX*HI*HJ**2*(HK+HI*SUM1~HJ*SUM2 )-HALFCF*( HI+HFI ) 

1  *< HI-HFI )*DSWDX/THETA 

D ( 1 ) =DTHRDX 
D ( 2 ) =DH I DX 
RETURN 
END 


SUBROUTINE  EMUAIR  <H,P,EMU) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  - -  1/30/71  -  W.J.GLOWACKI 

C 

C 

C  THIS  SUBROUTINE  EVALUATES  THE  VISCOSITY  OF  AIR  USING 

C  FITS  TO  THE  DATA  OF  HANSEN  ( NACA  TN  4150*  MAR  1968) 

C 

11  EMUSTR=(4.654E-08*SQRT  ( H ) ) / ( 1 . 0+48 . 47 /H ) 

12  IF(H-IOO.O) 13,13.15 

13  FMU=EMUSTR 

14  GO  TO  100 

15  IF(H-600.0)  16,16,19 

16  X=ALOG10 (H) -2.0 

17  EMU=EMUSTR*< 1.0-0. 066*X*X ) 

18  GO  TO  100 

19  Y=ALOG10(P/2116.0)+4.0 

20  Y2=Y*Y 

21  Y3= Y2*Y 

22  IF(H- 2500.0)23, 23, 32 

23  X=ALOG10(H)-2. 77815 

24  X2=X*X 

25  X3=X2*X 

26  A  =  2. 934 5-0.0 18 506* Y3  +  0. 341 56*Y2-1 .8930* Y 

27  B  =  -3. 0  797  +  0. 01 248 9* Y 3-0. 2 4747* Y2+1 .5311* Y 

28  C= 0.268 3 3-0. 00 16 71 3* Y3+0 . 0 32 4 5 5*Y2-0 . 19542*Y 

29  D  =  0 • 96  0 

30  GO  TO  99 

32  IF ( H-5000.0 ) 33 ,33  ,41 

33  X=ALOG10(H) -3. 39794 

34  X2=X*X 

35  X3=X2*X 

36  A=-15.510-0.31453*Y3+2.9963*Y2-4.3024*Y 

37  B=6. 172 2+0. 135 13* Y3-1 .3980* Y2+2. 6807* Y 

38  C=-0. 59 538-0. 010095*Y3+0. 13061 *Y2-0. 38801* Y 

39  0=0.642-0. 0005832 5*Y3+0. 0058 7 5*Y2+0 . 0 1 7082 *Y 

40  GO  TO  99 

41  IF1H-15000. 0)44,44*42 

42  WR I TE( 6  ,43  ) 

43  FORMAT ( 5  5H  H  EXCEEDS  15000  BTU/LB ,  SUBROUTINE  EMUAIR  IS  INVALID) 

44  X=ALOG10(H) -3.69897 

45  X2=X*X 

46  X 3  =  X 2*X 

47  A= 0.699 54+0. 02 1441 2* Y 3-0. 2172 l*Y2+0. 65 378* Y 

48  B=-0. 12 93 1-0. 005^775* Y3+0. 057 196* Y2-0.25586*Y 

49  C=-0. 42031-0. 001909*Y3+0.024503*Y2-0.0l3625*Y 

50  D=0.599+4. 175E-05*Y3+2. 500E-04*Y2+2. 5833E-02*Y 
99  EMU=EMUSTR*(A*X3+B*X2+C*X+D) 

100  RETURN 
END 
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SUBROUTINE  EMUN2  ( T E , ROE  * EMUE > 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  EVALUATES  THE  VISCOSITY  OF  NITROGEN  USING  THE 
BREBACH-THODOS  REDUCED-STATE  CORRELATION  FOR  DIATOMIC  GASES 
REF.  INDUSTRIAL  AND  ENGINEERING  CHEMISTRY  JULY  1968 

DATA  J/0/ 

TC  =  TEMPERATURE  AT  CRITICAL  POINT  (DEG. KELVIN) 

EMUTC  =  VISCOSITY  AT  1  ATM.  AND  TC  ( CENT  I PO I SES ) 

AO  =  CONVERSION  FROM  CENTIPOISES  TO  LBF  SEC/FT**2 
A 1  =  CONVERSION  FROM  G/CC  TO  SLUG/FT**3 

EVALUATE  CONSTANTS  AT  FIRST  ENTRY  OF  ROUTINE  ONLY 
IF  ( J.GT.C)  GO  TO  20 
TC  AND  EMUTC  FOR  NITROGEN 
10  TC=126.2 

EMUTC=8.65E-03 
A0=2.08854E-05 
Al=l .9403? 

A2=A0*EMUTC/TC**.979 
A3=1.7884*A0#EMUTC/SQRT ( TC) 

A4=0.7884*TC 

A5=1.196*A0*EMUTC/TC*#.659 

A6=.01198*A0/A1 

A7=.064! 0*A0/A1/A1 

T 1  =  TC 

T2=3.5#TC 

J  =  1 


COMPUTATION  OF  VISCOSITY  WITHOUT  DENSITY  CORRECTION 

BREBACH-THODOS  CORRELATON  FOR  T/TC  LESS  THAN  1 
20  IF  (TE.GT.T1)  GO  TO  30 
EMUE=A2*TE*#.979 
GO  TO  50 

SUTHERLAND-TYPE  LINK  FOR  T/TC  BETWEEN  1  AND  3.5 
30  IF  (TE.GT.T2)  GO  TO  40 

EMUE= A3*SQRT ( TE ) / ( 1 .+A4/TE ) 

GO  TO  50 

BREBACH-THODOS  CORRELATON  FOR  T/TC  GREATER  THAN  3.5 
40  EMUE=A5*TE**.659 


COMPUTATION  OF  DENSITY  CORRECTION  TO  VISCOSITY 
(NON-NEGATIVE  FIT  TO  BREBACH-THODOS  NITROGEN  DATA) 
50  TR=TE/TC 
D1=EMUE/A0 
D2  =  D1 /EMUT  C 
DMUE=0 . 

IF  ( ROE .LE. .0194  )  GO  TO  60 
DMUE=ROE*( A6+A7*ROE ) 

EMUE=EMUE+DMUE 
60  D3  =  ROE / A 1 
D4  =  DMU  E / AO 
RETURN 
END 
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FUNCTION  ENTHLP(U»X) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

C  THIS  SUBROUTINE  EVALUATES  H  FROM  LOGlO  P  AND  LOGIO  RHO 

C  AND  DIFFERS  FROM  ENTHLP  USED  IN  AIETP  DATED  7/01/65 

C  ONLY  IN  THAT  THE  ERROR  CORRECTIONS  ARE  NEGLECTED 

C 

AAl=2.4876733E-2+U* ( -8 . 2 5 53 89 0 E-3+U* ( 9 . 3 53 1 96 8 E-4+U* < -6 . 9 80 2 1 46 E-4 
l+U*(-7.0256556E-5+U*( 1 . 1 4492 0 8 E-4+U* ( 3 . 5 36 2 594E-5+U* ( 3 . 47 41 999 E-6 
2+8.9880020E-8*U)  )))))) 

AA2=3.4894155+U*( -2 . 9 1 9 3 9 1 1 E-2 +U* < 1 . 048 1 798 E-2+U* ( -2 . 671 473 3 E-4 
l+U*(-2.6951449E-4+U*(-2.31511 1 9 E-4+U* < -6 . 27498 93 E-5+U* (- 
27.5707012E-6-3.6330450E-7#U)  )))))) 

AO=.36547857+U*{ 1 . 0 1 8 92 1 0+U* ( -5 . 4844838 E-3+U* ( -2 . 7426243 E~2 +U* ( 
l-7.7107391E-3+U*< 2 . 42 0290 1 E-3 +U* < 1 . 425 6 5 73 E-3+U* < 2 . 2860249E-4+ 

21 .2352353E-5*U  ))))))) 

Al=-23. 014062+U* ( ,46035662+U* ( 3 . 004463 7E-2+U* ( 2 . 8 92 1 6 5 8 E-2+U* ( 

1- 4. 43371 02 E-3+U* ( 1 . 89 742 29E-4+U* ( 2 . 8 9 2 1472 E-3+U* ( 8 . 1 8 5 3 1 24E-4+ 
26.2730577E-5*U) )))))) 

A2=.71774971+U* ( -6 . 99 1 52 46E-2 +U* ( -4 . 767762 8 E-2+U* ( -7 . 2 32 8 36 7 E-3 
1+U*< 4.68771 00E-3+U* ( 2 . 1 1 95 796 E-3 +U* ( 2 . 6545 77 1 E-4+U* ( -4 . 641 6 05 1 E-7 

2- 1.4305115E-6+U) )))))) 

CO=.94993342+U* ( 1 . 0 1 8 73 00+U* ( - 2 . 78 3 1 274E-2 +U* ( -3 . 396820 8 E-2+U* ( 
1-7. 5758 964E-3+UD  4 . 33 77 1 4 5E-3 +U* ( 2 . 3 5 8 5 5 2 3 E-3+U* (  3 . 9680 5 2 9E -4+ 
22.2767082E-5+U) )))))) 

Cl=-56.5+44.6/ ( l.+EXP  (-U-1.751) 

C2=4.2598137+U*(-1.8016141+U*(  . 2 849 3 8 0 5+U* ( 2 . 5 5 240 3 4E-2+U* ( 

1- 2.806944CE-4+U*(-1.8133469E-3+U*(-9.7441476E-4+U*( - 1 . 8 3 565 1 1 E-4 

2- 1. 1383541E-5*U) )))))) 

EO=l. 374+1. 046*U 

El  =  -20.616464  +  U* ( . 99527941  +  U* (  .  36 3 1 8 074  +  U* ( 7 . 4026 95 4E-2  +  U# ( 

1-9. 5602994E-3+U* ( -6 . 868 3 405 E-3+U* ( -8 . 0 744407 E-4+U* ( 2 . 9434846E-5+ 
27.0238870E-6*U) )))))) 

E2=-.22188469+U*(-5.5351400E-2+U*(-8.3290583E-3+U*(-6.5111980E-4 
1+U*( 1 .3836166E-3+U* ( 8 . 2 097 1 36 E-4+U* ( 1 . 90 3 06 99 E-4+U* < 2 . 07022 3 9E- 5+ 
28.9265052E-7  +  U)  )))))) 

E3=.48990558+U* ( . 30389535+U+( 7 . 2 5 748 1 7E-2 +U* ( -6 . 9 7486 1 5E-3+U* ( 

1- 6.4693200E-3+U*( -5 . 8 78 87 07 E-4+U* ( 7 . 7 2 3262 9 E-5+U* (5.3091414E-6- 
26.6747741E-7+U ))))))) 

GO=l. 835+1 .043*U 

Gl=-11.591952  +  U*( 1 .8211028+U* ( 1 . 868 2006 E-2 +U* ( 3 . 8 6275 52 E-2  +  U*  ( 
11.7205973E-2+U*(-3.7323972E-3+U*(-2.2956386E-3+U*(-3.076423'9E-4- 
21.2352353E-5  +  U)  )))))) 

G2  =  1.348 070 5  +  U*  (  .  2 1  904060  +  U*  <  -6. 7  5  5242  9E-3+.U*  (  1 . 6  5 6 34 52 E-2  +  U*  ( 

19. 1140456E-3+U*(-2.3675397E-3+U*(-1.7376722E-3+U# (-2.9616734E-4- 
21.6477373E-5+U)  )))))) 

G3=-2. 26283 66+U*( -1 .8824979+U* (-. 28388494+U* ( 1 . 2 9 0 5 52 5 E-2+U* ( 
11.4367226E-2+U*(6.8502795E-6+U*(-9.3289221E-4+U*( - 1 . 7 5 5 1 8 92 E-4 

2- 1.0293628E-5+U ))))))) 

HO=2. 376+1 .046*U 
Hl=-28 .86+. 507+U 

H2=.95+.38/(  li  +  EXP  (2.5+IU+4.2) ) ) 

Z1=E1*( X-EO) 

IF1Z1-10. ) 11 .10,10 

10  FX 1  =0. 0 
GO  TO  21 

11  IF(Zl+8. (12,12.13 

12  FX1=E2*X+E3 
GO  TO  20 

13  FXl=(E2*X+E3)/( l.+EXP  (Zl)) 

20  Z2  =  G 1* ( X-GO ) 

IF(Z2-10. >22,21,21 

21  FX2  =  0 . 0 
GO  TO  31 
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22  IFIZ2+8.  >23,23*24 

23  FX2=G2*X+G3 
GO  TO  30 

24  FX2= (G2*X+G3 ) / ( 1 .+EXP  (Z2)) 

30  Z3  =H 1* ( X-HO ) 

IF(Z3-10. >32,31,31 

31  FX3=0.0 
GO  TO  40 

32  IFIZ3+8. >33,33,34 

33  FX3=H2*(X-HO) 

GO  TO  40 

34  I F { ABS  (Z3)-. 01)35  *35, 36 

35  FX  3  =  -H2 /H 1 
GO  TO  40 

36  FX3=H2* < X-HO) / ( 1 .-EXP  ( Z3 )  ) 

40  XX=X-FX1-FX2-FX3 

50  Z4=A1*< XX-AO) 

IF(74-10. >52,51,51 

51  F 1 =0 . 0 
GO  TO  61 

52  IF(Z4+8. )53,53*54 

53  F1=A2*< XX-AO) 

GO  TO  60 

54  I F ( ABS  (Z4)-. 01 >55,55,56 

55  F 1  =  -A 2/ A 1 
GO  TO  61 

56  F1=A2*( XX-AO) /( 1 .-EXP  (Z4)) 

60  Z5=C1*( XX-CO) 

IF1Z5-10. >62,61,61 

61  F2  =  0 . 0 
GO  TO  70 

62  I F< Z5+8  . ) 63 ,63 ,64 

63  F2=C2*< XX-CO) 

GO  TO  70 

64  I F ( ABS  (Z5)~. 01 >65,65,66 

65  F2=-C2/C1 
GO  TO  70 

66  F2=C2*(XX-CO)/(l.-EXP  (Z5)> 

70  HRHOVP  =  AA 1*XX+AA2+F1+F2 

75  ENTHLP=33.687746*HRH0VP*10.** ( X-U) 
RETURN 
END 


•  SUBROUTINE  GAUSS 
C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

C 

C  THIS  SUBROUTINE  EVALUATES  DELSTR/THETA  AND  DELTA/THETA  BY 

C  16  POINT  GAUSSIAN  INTEGRATION 

C 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT,XT,HO, 

1  DXDX»DRDX,DEMPX, DRHODX  ,DPEDX,DUEDX,DHEDX* 

2  CF,COEFF,DELSTR ,DELTA,DHIDX ,DH1 ,DH2 ,DRWDX, 

3  DSVDEL,DSVTH,DSWDX,DTHDX,EMUE,GAS,HAD,HFI  , 

4  HI ,HJ,HK,HW, IP, I  PI ,LAWCF,LM,M,N,PAGE,PR,Q> 

5  RE , RECVRY  »  RETH  »REY  ANL , RHOW  » ROE , RW ,ST , SUM1 , 

6  SUM2,TE, THETA, THVDEL,TW,TX(8)  ,  XEND , XPUNCH , 

7  XTRA1»XTRA2»XTRA3,ZE 
DIMENSION  C ( 16 ) ,Z ( 16) 

C 
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C 


C 


COEFFICIENTS  FOR  16  POINT  GAUSSIAN  INTEGRATION 
DATA  C/  .01357623*  .031126762*  .047579256*  .062314485* 

.08457826,  .091301710,  .094725305,  .094725305, 

.08457826,  .074797995,  .062314485,  .047579256, 


1 

2 

3 

4 

5 

6 


.01357623/ 


,Z/  .005299535,  .027712490, 


.12229780,  .191061880,  .270991610,  .359198220, 
.54750626,  .640801780,  .729008390,  .808938120, 
.93281560,  .972287510,  .994700460/ 


.074797995, 

.091301710, 

.031126762, 

.067184400, 

.452493740, 

.877702200, 


SUM1=0. 0 
SUM2=0 . 0 
SUM3=0 . 

SUM4=0 . 

R0  =  R0E 

DO  10  J=1 , 1 6 
1  =  1 7- J 

UVUE  =  Z (  I  )**HK 
H=HW+UVUE*( DH1-UVUE*DH2  ) 

CALL  RHOTZ  ( GAS , P E , H , RO , T , Z EE ) 
TERM=(ROE/RO-l» )*C(  I  ) 
TERM3=UVUE*Z ( I ) * ( ROE/RO ) *C ( I > 
SUM1=SUM1+TERM 
SUM2=SUM2+TERM*Z(  I  ) 
SUM3=SUM3+TERM3 
SUM4=SUM4+TERM3*UVUE 
10  CONTINUE 
A=SUMl + 1 . 

XTR A2= ( 1.-SUM3/A) /A 
XTRA3= ( 1 .-SUM4/A) /A 
COEFF= ( HI+1 . ) *H I / { HI -1 . ) 
DSVTH=COEFF*SUMl+HI 
DELVTH=COEFF*( SUM  1  +  1 . ) 

THVDEL= l./DELVTH 
DSVDEL=DSVTH*THVDEL 
RETURN 
END 


SUBROUTINE  HWALL  ( GAS , P , T , H ,RH0 ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  COMPUTES  HW  FOR  A  GIVEN  PE  AND  TW 

DATA  C/2.3025851/,  RHO0/2.423516E-03/ ,  AIR/6HAIR 

IF  (GAS. NE. AIR)  GO  TO  10 

AIR  CALCULATION 

FOR  AIR,  WE  ASSUME  PRESSURE  IS  LOW  ENOUGH  SO  HW  =  HW(T)  ONLY 
DT  =  T-499 • 17 

H=.123406*T*(3.48+.00033175*DT/< 1 . -EXP < - . 012287*DT > ) ) 

GO  TO  60 

NITROGEN  CALCULATION 
10  Z= 1  • 

RHOZ=P/ (3196. 8*T) 

IF  (RHOZ.LE.RHOO)  GO  TO  40 
DO  20  1=1,50 
X  =  ALOG1 0 ( RHO /RHOO ) 

Y=.00482*X*EXP(X) 

IF  (X.LE.O.)  Y=0 . 

Z  =  1 . +Y*  X 
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F=RHO*Z-RHOZ 
F1=Z+Y* (X+2. ) /C 
RH0=RH0-F/F1 
Q=ABS ( -F/Fl /RHO) 

IF  (Q. LE.. 000001 )  GO  TO  50 
20  CONTINUE 

WRITE  (6.30) 

30  FORMAT  (28H0  NO  CONVERGENCE  IN  HWALL  ) 

40  RHO=RHOZ 

50  H=.  12772*T*( Z  +  2.5  ) 

IF  (T.LE.500. )  GO  TO  60 
B=EXP(-3390./T) 

H=H+432.971#B/ ( 1 . — B ) 

60  RETURN 
END 


SUBROUTINE  INTERP  ( J1 .JMAX.X.Y .XO.YO.JJ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  PERFORMS  THREE  POINT  INTERPOLATION  IN  A  TABLE 
OF  X  AND  Y  TO  FIND  YO  FOR  A  GIVEN  XO.  THE  SEARCH  IS  IN  ONE 
DIRECTION  ONLY.  BEGINNING  AT  J1  AND  ENDING  AT  END  OF  TABLE. 

DIMENSION  X< JMAX)  .Y< JMAX ) 

DATA  K/l/ 

GO  TO  ( 10.20)  .  K 
10  K  =  2 

SI GN=+1  . 

IF  (X(  Jl).GT.X(Jl-l) )  GO  TO  20 
SIGN=-1. 

20  SX0=SIGN*X0 

DO  30  J  =  J 1 » JMAX 
SXAV=SIGN*(X( J-l )+X( J) ) /2. 

IF  (SXO.LT.SXAV)  GO  TO  40 
30  CONTINUE 
40  DX01=X0~X( J-2) 

DX2 1 =X ( J-l ) -X ( J-2 ) 

DX31  =  X( J!-X ( J-2  ) 

DY2 1 =Y ( J-l } - Y ( J-2  ) 

DY3 1  =  Y ( J ) -Y ( J-2  ) 

B=( DY31/DX31-DY21/DX21 ) / ( DX31-DX21  ) 

A=DY21/DX21— B*DX21 
YO=Y( J-2)+DX01*( A+B*DX01  ) 

J J= J-2 
RETURN 
END 


SUBROUTINE  LABEL1  (GAS) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  — - -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  READS  A  MAXIMUM  OF  50  CARDS  CONTAINING  HOLLORITH 
INFORMATION  WHICH  IS  PRINTED  AT  BEGINNING  OF  OUTPUT.  FOR  CASES 
AFTER  THE  FIRST  CASE,  THIS  ROUTINE  PERMITS  LINES  1  THROUGH  LIN 
PLUS  UP  TO  10  ADD  I T I T I ONAL  LINES  (ARBITRARILY  DISTRIBUTED  AMONG 


C-21 


nnooooo  r>  ooo  o  o  o 


NOLTR  71-6 


C  THE  FIFTY  ALLOWED)  TO  BE  REPLACED  BY  NEW  LINES  INPUT  FROM  CARDS. 

DIMENSION  LREPLC( 10 ) ,TITLE( 600  ) 

DATA  AIR/6HAIR 

READ  IDENTIFYING  AND  DESCRIPTIVE  TITLING  (50  LINES  MAXIMUM) 

10  READ  (5,20)  LIN, LOUT, (LREPLC(L),L=1*10) 

20  FORMAT  (1215) 

I  I N  =  12*L I N 
I  OUT  = 1 2*LOUT 
IF  (LIN.EQ.O)  GO  TO  40 
READ  (5,30)  (TITLE!  I)  ,1  =  1, IIN) 

30  FORMAT  (12A6) 

40  DO  50  L  =  1 ,10 

IF  (LREPLC(L) .EQ.O)  GO  TO  50 
I12=12*LREPLC(L) 

I  1=112-11 

READ  (5,30)  (TITLE!  I ) ,1  =  11,112) 

50  CONTINUE 


WRITE  IDENTIFYING  AND  DESCRIPTIVE  TITLING  (50  LINES  MAXIMUM) 

60  IF  (GAS. EQ. AIR)  GO  TO  80 
WRITE  (6,70) 

WRITE  (16,701 

70  FORMAT  (1H1  27X , 45HN I TROGEN  TURBULENT  BOUNDARY  LAYER  CALCULATION 

*  /28X  »45H - /////) 

GO  TO  100 

80  WRITE  (6,90) 

WRITE  (16,90) 

90  FORMAT  (1H1  27X,40HAIR  TURBULENT  BOUNDARY  LAYER  CALCULATION 

*  /28X  »40H -  /////) 

100  IF  (IOUT.EO.O)  GO  TO  130 

WRITE  (6,110)  ( T I TLE ( I  )  ,1  =  1,1  OUT) 

WRITE  (16*110)  (TITLE! I ) ,1=1, IOUT) 

110  FORMAT  (1H  27X*12A6) 

WRITE  (6,120) 

WRITE  (16,120) 

120  FORMAT  ( 1H0  //  43X,43H(THE  NEXT  PAGE  IS  A  DUPLICATE  OF  THIS  PAGE)) 
130  IF  (GAS. EQ. AIR)  GO  TO  140 
WRITE  (6,70) 

WRITE  (16,70) 

GO  TO  150 
140  WRITE  (6,90) 

WRITE  (16,90) 

150  IF  (IOUT.EQ.O)  GO  TO  200 

WRITE  (6,110)  (TITLE!  I )  ,1=1, I  OUT) 

WRITE  (16,110)  ( T I TLE (  I  ) , I  =  1 , I  OUT ) 

200  RETURN 
END 


SUBROUTINE  LABEL2  ( XO , RO , OMEG AO , NN , J ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  WRITES  AND  IDENTIFIES  FOR  RECORD  PURPOSES  THE 
DATA  INPUT  IN  THE  MAIN  ROUTINE  AND  THE  TEMPERATURE  DISTRIBUTION 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT ,XEXIT  ,XT  *H0» 

1  DXDX,DRDX,DEMDX,DRHODX ,DPEDX,DUEDX ,DHEDX, 
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2  CF,COEFF»DELSTR, DELTA, DHIDX,DH1,DH2»DRWDX» 

3  dsvdel»dsvth»dswdx,dthdx,emue»gas»had,hfi, 

4  HI ,HJ,HK,HW, IP, IP1,LAWCF,LM,M,N,PAGE,PR.Q, 

5  RE,RECVRY,RETH,REYANL, RHOW ,ROE,RW,ST ,SUM1 , 

6  SUM2»TE,THETA,THVDEL,TW,TX(8) , X END , XPUNCH , 

7  XTRA1.XTRA2.XTRA3.ZE 
DATA  AIR/6HAIR  / 

C 

C 

IF  (GAS. EQ. AIR)  GO  TO  30 

10  WRITE  (  6,20)  NN, THETA, REXIT, XPUNCH, X0.R0.OMEGA0, HI ,LM, IP, J.TX 
WRITE  (16,20)  NN, THETA, REXIT, XPUNCH, XO » RO » OMEGAO , HI  ,LM, IP, J.TX 
20  FORMAT  ( 1 HO /41 HON  I TROGEN  TURBULENT  BOUNDARY  LAYER  NUMBER  15// 

*  8H0THETA  =  1PE14 . 7 , 8X , 7HREX I T  =  1PE14.7 ,9X  .9HXPUNCH  =  1PE14.7/ 

*  8H0  XO  =  1PE14.7,10X»6HR0  =  1PE14 . 7 , 9X , 9HOMEGA0  =  1PE14.7/ 

*  8H0  HI  =  1PE14.7,11X,4HLM  =  I2,11X,4HIP  =  I2.10X.3HJ  =,I2// 

*  2  5H0TEMPERATURE  DISTRIBUTION  /7X  ,  2HX* „ 1  OX  ,  2HT*  1  OX , 2HX 1 ^ 1 OX , 2HT 1 

*  1  OX , 2HX2  1  OX , 2HT 2  10X.2HXE  10X.2HTE/  4 ( OPF 1 2 . 6 , F 1 2 . 4 )  //  ) 

GO  TO  50 

30  WRITE  (  6,40)  NN , THET A , REX  I T , XPUNCH , XO , RO , OMEGAO ♦ H I , LM , I P , J  ,  TX 
WRITE  (16,40)  NN, THETA, REXIT, XPUNCH.XO.RO, OMEGAO, HI  ,LM, IP, J.TX 
40  FORMAT  (1H0/36H0AIR  TURBULENT  BOUNDARY  LAYER  NUMBER  15// 

*  8H0THETA  =  1 PEI 4. 7 , 8X , 7HREX I T  =  IP E 14 . 7 , 9X , 9HXPUNCH  =  1PE14.7/ 

*  8H0  XO  =  1PE14.7»10X,6HR0  =  1PE14.7  »9X,9H OMEGAO  =  1PE14.7/ 

*  8H0  HI  =  1PE14.7,11X,4HLM  =  I2.11X.4HIP  =  I2,10X,3HJ  =  12// 

*  25H0TEMPERATURE  DISTRIBUTION  /7X.2HX*  10X.2HT*  10X.2HX1  10X.2HT1 

*  10X.2HX2  10X.2HT2  10X.2HXE  10X.2HTE/  4 ( OPF 1 2 . 6  ,  F 1 2 . 4 )  //  ) 

C 

50  GO  TO  (60,80)  ,  LAWCF 
60  WRITE  (  6,70) 

WRITE  (16,70) 

70  FORMAT  ( 46H0LAWCF  =  1  ( SPALD I NG-CH I  SKIN  FRICTION  LAW)  ) 

GO  TO  100 
80  WRITE  (  6,90) 

WRITE  (16,90) 

90  FORMAT  ( 80H0LAWCF  =  2  ( LUDW I EG-T I LLM ANN  SKIN  FRICTION  LAW 

*  WITH  REFERENCE  ENTHALPY)  ) 

C 

100  RETURN 
END 


SUBROUTINE  LABEL3  ( LM ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  WRITES  A  LIST  IDENTIFYING  THE  HEADINGS  USED 
FOR  THE  TABULATED  RESULTS  OF  THE  BOUNDARY  LAYER  CALCULATION 

WRITE  (6,10) 

10  FORMAT  (1H1////44X.41HIDENTIFICATI0N  OF  HEADINGS  AND  UNITS  USED 

*  /44X.41H - /// 

*//38X,55HD — /DX  DENOTES  THE  DERIVATIVE  OF  —  WITH  RESPECT  TO  X  // 
*//38X,57HCDLSTR  AXISYMMETRIC  B.L.MASS  CONSERVATION  TH I CKN ESS ( F T  ) 
*/ / 38X , 57HCF  SKIN  FRICTION  COEFFICIENT 

*//38X,57HCUTOFF  OPTIMUM  CUTOFF  POINT  OCCURS  WHEN  CUTOFF  =  0 
*// 38X , 57HD — /DX  DERIVATIVE  OF  —  WITH  RESPECT  TO  X 
*//38X,57HDELSTR  2-D  BOUNDARY  LAYER  DISPLACEMENT  THICKNESS  (FT) 
*//38X,57HDELTA  BOUNDARY  LAYER  THICKNESS  (FT) 

*//38X,57HDLINCH  BOUNDARY  LAYER  THICKNESS  (INCH) 

* / / 38X , 5 7HDLSTRR  AXISYMMETRIC  BOUNDARY  LAYER  THICKNESS  (FT) 
*//38X,57HDSVDEL  DELSTR  /  DELTA 

'  *//38X,57HDSVTH  DELSTR  /  THETA 
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*/'/38X*57HEM  MACH  NUMBER  AT  EDGE  OF  BOUNDARY  LAYER 

*//38X»56HEMUE  VISCOSITY  AT  EDGE  OF  B.L.  ( LBF  SEC/FT2  )  ) 

WRITE  (6,20) 

20  FORMAT  ( 

*  / 38X , 5 7HERROR ( 1 )  TRUNCATION  ERROR  IN  INTEGRATION  OF  DTHETA/DX 
*//38X,5 7HERROR ( 2 )  TRUNCATION  ERROR  IN  INTEGRATION  OF  DHI/DX 

* / / 38X , 5 7HHAD  ADIABATIC  WALL  ENTHALPY  ( BTU/LBM ) 

*//38X,5 7HHE  ENTHALPY  AT  EDGE  OF  BOUNDARY  LAYER  (BTU/LBM) 

* / / 38 X  »  5  7HHF I  SHAPE  PARAMETER  IN  TRANSFORMED  PLANE  FOR  DP/DX=0 

*//38X,57HHI  SHAPE  PARAMETER  IN  TRANSFORMED  PLANE 

*//38X  ,5  7HHW  WALL  ENTHALPY  (BTU/LBM) 

*/ / 38X  *  5 7HPE  PRESSURE  AT  EDGE  OF  BOUNDARY  LAYER  ( LBF/FT2 ) 

*//38X,5 7HQ  HEAT  TRANSFER  RATE  AT  NOZZLE  WALL  ( BTU/FT2*SEC ) 

*/lHl/// 

*//38X,57HR  RADIAL  DISTANCE  FROM  AXIS  TO  EDGE  OF  B.L.  (FT) 

*//38X,57HRE  REYNOLDS  NUMBER  PER  FOOT  AT  EDGE  OF  B.L. 

*//38X,56HRETH  REYNOLDS  NUMBER  BASED  ON  THETA  ) 

WRITE  (6,30) 

30  FORMAT  ( 

*  / 38X , 5  7HRHOE  DENSITY  AT  EDGE  OF  BOUNDARY  LAYER  ( SLUG/FT3 ) 

#//38X,57HRINCH  RADIAL  DISTANCE  FROM  AXIS  TO  EDGE  OF  B.L.  (INCH) 

*/ / 3 8X  *  5  7HROE  APPROXIMATE  DENSITY  AT  EDGE  OF  B.L.  ( SLUG/FT  3 ) 

*//38X,57HRW  RADIAL  DISTANCE  FROM  AXIS  TO  NOZZLE  WALL  (FT) 

*//38X,57HRWINCH  RADIAL  DISTANCE  FROM  AXIS  TO  NOZZLE  WALL  (INCH) 

*//38X,57HST  STANTON  NUMBER  CALCULATED  USING  REYNOLDS  ANALOGY 

*/ / 38 X  »  5  7HSW  DISTANCE  ALONG  NOZZLE  WALL  CONTOUR  (FT) 

*// 38X  ,5  7HTE  APPROXIMATE  TEMPERATURE  AT  EDGE  OF  B.L.  (DEG  K) 

*//38X,57HTHETA  2-D  BOUNDARY  LAYER  MOMENTUM  THICKNESS  (FT) 

*//38X,57HTHETAR  AXISYMMETRIC  B.L.  MOMENTUM  THICKNESS  (FT) 

*//38X,57HTW  WALL  TEMPERATURE  (DEG  K) 

* / / 38X , 5 7HUE  VELOCITY  AT  EDGE  OF  BOUNDARY  LAYER  (FT/SEC) 

*// 38X  » 5 7HX  AXIAL  DISTANCE  FROM  NOZZLE  THROAT  (FT) 

*//38X,57HXINCH  AXIAL  DISTANCE  FROM  NOZZLE  THROAT  (INCH) 

*//38X,57HZE  APPROXIMATE  COMPRESSIBILITY  AT  EDGE  OF  B.L. 

*////38X ,53HD — /D X  DENOTES  THE  DERIVATIVE  OF  —  WITH  RESPECT  TO  X  ) 
C 

IF  (LM.EQ.O)  GO  TO  40 
WRITE  (16,10) 

WRITE  (16,20) 

WRITE  (16,30) 

C 

40  RETURN 
END 


SUBROUTINE  OUT  ( X , Z , D , ERROR , K , L , H ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  WRITES  THE  OUTPUT  OF  THE  BOUNDARY 
LAYER  CALCULATION  IN  THE  NOZZLE  ON  THE  SYSTEM  TAPE 
AND,  IF  LM  IS  GREATER  THAN  ZERO,  ON  TAPE  16  ALSO. 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT,XT,HO, 

1  DXDX,DRDX,DEMDX,DRHODX , DPEDX , DUEDX , DHEDX , 

2  CF , COEFF ,DELSTR, DELTA, DHIDX ,DH1 »DH2 ,DRWDX, 

3  DSVDEL,DSVTH’,DSWDX,DTHDX,EMUE,GAS,HAD,HFI  , 

4  HI ,HJ,HK,HW, IP , I  PI ,LAWCF,LM,M,N , PAGE, PR, Q, 

5  RE,RECVRY,RETH,REYANL,RHOW,ROE,RW,ST,SUMl, 

6  SUM 2 ,TE,THETA,THVDEL,TW,TX ( 8  )  , X END , XPUNCH , 

7  XTRA1 ,XTRA2,XTRA3,ZE 

8  ,CDLSTR,DLSTRR»THETAR,DTHRDX 
DIMENSION  Z(5) »D(5) ,ERROR(30) 
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DATA  XPREV/-100./ 

IF  (X.EQ.XPREV)  GO  TO  70 
XPRE V=X 
IP1=IP1+IP 
IF  (IP1.EQ.1)  1=0 

IF  (I.GT.O)  GO  TO  30 

1=7 

PAGE=PAGE+1  . 


WRITE  (6,10)  PAGE 
10  FORMAT  ( 1H1  111X,4HPAGE  F4.0//> 


WRITE 

(  6,20  ) 

20  FORMAT 

( 12OH0 

X 

RW 

R 

E 

*M 

RHOE 

PE 

UE 

hE 

*/ 1 20H 

DSW/DX 

DRW/DX 

DR/DX 

DEM/DX 

*  DRHOE/DX 

DPE/DX 

DUE/DX 

DHE/DX 

*/120H 

TW 

HW 

HAD 

TE 

* 

ZE 

ROE 

ERROR! 1 ) 

ERROR! 2) 

*/120H 

CUTOFF 

DSVDEL 

DSVTH 

XINCH 

*  RW INCH 

RI  NCH 

DLSTRR 

CDLSTR 

*/120H 

RE 

RETH 

ST 

EMUE 

*  DH I /DX 

DTHETAR/DX 

DELSTR 

DELTA 

*/ 120H 

POINT 

CF 

Q 

HFI 

* 

HI 

THETAR 

THETA 

DLINCH 

) 

IF  (LM.EQ.O)  GO  TO  30 
WRITE  (16,10)  PAGE 
WRITE  (16,20) 

C 

30  1=1-1 

CUTOFF=12.*(RW-DELTA-REXIT*<X-XT)/(XEXIT-XT)  ) 

XINCH=12.*X 

RWINCH=12.*RW 

RINCH=12.*R 

DSINCH=12.*DELSTR 

DLINCH=12.*DELTA 

WRITE  (6,40)  X,RW,R,EM,RHOE,PE,UE*HE,DSWDX,DRWDX,DRDX,DEMDX, 

1  DRHODX ,DPEDX,DUEDX,DHEDX , TW , HW , HAD , TE , ZE , ROE , ERROR ( 1 )  , 

2  ERROR! 2 ) , CUTOFF , DSVDEL , DSVTH , X I NCH , RW I NCH , R I NCH , DLSTRR , CDLSTR , RE  , 

3  R ETH , ST ,EMUE,DHIDX,DTHRDX,DELSTR,DELTA,IP1,CF,G,HFI ,HI ,THETAR, 

4  THETA  *  DL I NCH 

40  FORMAT  ( 1H  / 5 ( 1 P8 El  5 . 6 / )  1 1 0 , 5X ,  1 P7 E 1 5 . 6 ) 

IF  (LM.EG.O)  GO  TO  50 

WRITE  (  16»40)X,RW,R»EM,RHOE,PE,UE*HE*DSWDX»DRWDX*DRDX»DEMDX  , 

1  DRHODX  » DP EDX, DUE DX  »DHEDX , TW, HW,HAD , TE , Z E , ROE , ERROR  (  1 )  , 

2  ERROR! 2) , CUTOFF, DSVDEL, DS V TH, X.  INCH, RW INCH, RI NCH, DLSTRR, CDLSTR, RE, 

3  R ETH, ST ,EMUE ,DHIDX ,DTHRDX ,DELSTR , DELTA ,IP1,CF»Q,HFI,HI  ,  THETAR , 

4  THETA, DLINCH 

50  IF  ( X.LT.XPUNCH )  GO  TO  70 

PUNCH  60»X,RW»R,XI NCH , RW I NCH 
60  FORMAT  ( OPF 1 0 . 5 , 1 P4E 1 5 . 7 ) 

70  RETURN 
END 


SUBROUTINE  RHOTZ  ( GAS , P , H , RHO , T , Z  ) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  - -  W.J.GLOWACKI 

C 

C  THIS  SUBROUTINE  COMPUTES  RHO,  T,  AND  Z  FOR  A  GIVEN  P  AND  H. 

C  FOR  NITROGEN,  THE  CURVE-FITS  OF  HARRIS  ARE  USED,  WHILE 

C  FOR  AIR,  THE  CURVE-FITS  OF  GRABAU  ARE  USED. 

C 
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Data  c/2.3025851/ ,  RHO0/2.423516E-03/,  AIR/6HA1R  / 


IF  (GAS. EQ. AIR)  GO  TO  60 

NITROGEN  CALCULATION 
PATM=P/21 16. 

HOVR=H/. 127716 
ROAMAG=RHC/RHOO 
DO  30  1=1.50 

X  =  ALOG 1 0 ( ROAM AG ) 

IF  ( X.GT.O.  )  GO  TO  10 
Y  =  0. 

Z  =  l. 

T  =  HOVR / 3.5 

IF  ( T.GT.500. )  GO  TO  20 
ROAMAG=273. 16*PATM/T 
GO  TO  50 

10  Y=.00482*X*EXP(X) 

Z=1 .+Y*X 

20  T=273.16*PATM/(Z*ROAMAG) 

E0=3390. /T 
E  1  =  0  . 

IF  (T.GT.500.)  El=EXP(-3390./T! 

IF  (E0.LE..01)  E2  =  T/  ( 1 . -E0* ( .5-E0/6.  )  ) 

IF  (E0.GT..01)  E2  =  3390./(  1.-E1  ) 

F  =  HOVR- ( Z  +  2 .5 ) *T-E1*E2 

F 1  =  T* ( Z+( 2. 5+El*< E2/T )**2 )*< 1 .+Y* ( X  +  2 . ) /Z/C)  ) ZROAMAG 

ROAMAG=ROAMAG-F/Fl 

Q=ABS(-F/F1/R0AMAG) 

IF  (Q.LE.. 000001)  GO  TO  50 
30  CONTINUE 

WRITE  (6.40)  Q 

40  FORMAT  ( 28H0  NO  CONVERGENCE  IN  RHOTZ.  3HQ  =F9.7) 
50  RHO=ROAMAG*RHOO 
GO  TO  130 

AIR  CALCULATION 
60  X4=ALOG1C(RHO/0, 0025089296) 

X=ALOG10(P/2116,0) 

Y4=H-ENTHLP(X4»X) 

X2  =  X4 
Y2  =  Y4 

X4=0.95*X4 
X3  =  X4 

Y4=H-ENTHLP ( X4.X ) 

Y3  =  Y4 

DX=X3-X2 

DY= Y3- Y2 

X4=X3-Y3*DX/DY 

Y4=H-ENTHLP(X4»X) 

DO  100  1=1,20 

XI  =  X2 
Y 1  =  Y  2 
X2  =  X3 
Y2  =  Y3 
X3  =  X4 
Y3  =  Y4 
DXP=DX 
DYP=DY 
DX  =  X3- X  2 
D Y  =  Y3-Y  2 

DI V=DXP*DY-DX*DYP 
I F ( D I V . EQ . 0 .  )  GO  TO  70 
DEL INV= (DX+DXP)*DXP/DIV 

RAD  I C= ( DX  +  DY*DELINV)**2-4.*Y3^DX*DELINV 
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IFIRADIC.LT. 0.  )  60  TO  80 
XX=SQRT  (RADIO 
XL=X3+X2-DY*DELINV 
X4  = . 5* ( XL+XX) 

XL= • 5* ( XL-XX ) 

XBAR=X3-Y3*DX/DY 
DXR= ABS  ( X4-XBAR ) 

DXL=ABS  (XL-XBAR) 

IF(DXR.LE.DXL)  GO  TO  90 
X4  =  XL 
GO  TO  90 

70  X4=X3— Y3*DX/DY 
GO  TO  90 

80  X4=Y3-.25*(DX+DY*DELINV)**2/(DELINV*Dx) 

90  Y4=H-ENTHLP( X4,X) 

TEST=ABS  ( ( X4-X3 ) /X4 ) 

IF  (TEST. GT.. 0001  )  GO  TO  100 
IF  ( Y4.GE.-.001.AND.Y4.LE. .001 )  GO  TO  120 
100  CONTINUE 

WR I TE ( 6  » 1 10 ) 

110  FORMAT  ( 28H0  NO  CONVERGENCE  IN  RHOTZ,  6HTEST  =  F7.5) 
120  R=0.0025089296*EXP  ( X4/ 0 . 43429448 ) 

RHO  =  R 

Z=C0MP ( X4»X) 

T=3.2373132E-04*P/Z/RHO 

C 

130  RETURN 
END 


SUBROUTINE  RHOTZE 
C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  - 1/30/71 - W.J.GLOWACKI 

C 

C 

C  THIS  SUBROUTINE  EVALUATES  THE  DENSITY,  TEMPERATURE,  AND 

C  COMPRESSIBILITY  ( ROE ,  T  E  ,  ZE  )  AT  THE  EDGE  OF  THE  BOUNDARY 

C  LAYER.  FOR  NITROGEN,  THE  VALUES  ARE  ONLY  APPROXIMATE  IF  THE 

C  I SENTROP I C  EXPANSION  DATA  USED  IN  DETERMINING  THE  INVISCID 

C  CORE  WAS  OBTAINED  FROM  WOOLLEYS  PROGRAM.  IF  THE  DATA  WAS 

C  OBTAINED  FROM  THF.  EXPANSION  PROGRAMS  LABELED  AIETP  OR  NIETP, 

C  THEN  THE  VALUES  ARE  EXACT  AS  THE  SAME  EQUATIONS  ARE  USED. 

C 

COMMON  EX,R,EM»RHOE,PE,UE,HE,REXIT  »XEXIT  »XT  ,H0» 

1  DXDX ,DRDX ,DEMDX,DRHODX,DPEDX,DUEDX ,DHEDX  , 

2  CF,COEFF,DELSTR, DELTA, DHIDX,DH1,DH2,DRWDX, 

3  DSVDEL,DSVTH,DSWDX,DTHDX,EMUE,GAS,HAD,HFI  , 

4  HI ,HJ»HK*HW»IP»IP1,LAWCF ,LM,M,N,PAGE»PR,Q, 

5  RE,RECVRY,RETH,REYANL,RH0W,R0E,RW,ST,SUM1 , 

6  SUM 2  ,  TE , THETA, THVDEL ,TW  ,  T  X ( 8 ) , XEND  *  XPUNCH , 

7  XTRA1 ,XTRA2»XTRA3*ZE 

DATA  AIR/6HAIR  / , RHOO , P 0 / . 00 2 5 089 296 , 2 1 1 6 . / 

C 

C 

ROE=RHOE 

IF  (GAS. EC. AIR)  GO  TO  10 
C 

C  NITROGEN  CALCULATION 

CALL  RHOTZ  ( GAS  , P E , HE , ROE , T E , Z E ) 

GO  TO  20 
C 

C  AIR  CALCULATION 
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10  U=ALOG10(RHOE/RHO0l 
X  =  ALOG10( PE/PO  ) 

ZE  =  COMP ( U  j  X ) 

TE=3.2373132E-04*PE/ZE/RHOE 

20  RETURN 
END 


SUBROUTINE  SHAPE  ( X0*Y0,DYDX0*X»Y  » JMAX  *K0  » I ) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

C 

C  THIS  SUBROUTINE  DETERMINES  THE  SHAPE  OF  THE  SUBSONIC 

C  INLET  DEFINED  BY  QUARTIC  Y= ( ( ( A*X+B ) *X+C ) *X+D ) * X+E 

C  USING  R,  DR/DX  SPECIFIED  AT  BEGINNING  OF  INLET  <X=X0>* 

C  AND  USING  R,  DR/DX  (=0).  D2R/DX2  AT  NOZZLE  THROAT  (X=0> 

C  WITH  D2R/DX2  DETERMINED  BY  FITTING  A  QUADRATIC  (WITH  ZERO  SLOPE 

C  AT  THROAT)  TO  THROAT  AND  FIRST  CHARACTERISTIC  CORE  POINT 

C 

DIMENSION  X(JMAX) .Y(JMAX) 

C 

X03=X0**3 
YST AR  =  Y (  I  ) 

DETERMINE  CONSTANTS  FOR  QUARTIC  INLET  GEOMETRY 
E=YSTAR 
D  =  0. 

C=  (  Y  ( 1+1 )-YSTAR)/X( I +1 )  **2 
B= ( 4 « * ( Y0-YSTAR)-X0*(2.*C*X0  +  DYDX0)  )/X03 
A=(  (C*X0+DYDX0)-3.*<  YO-YSTAR) /X01/X03 

INITIAL  GUESS  FOR  EX  NEAR  THROAT  FROM  TAYLOR  SERIES  EXPANSION 
Y-Y*  =  ( X-X*) ( DY/DX )*  +  0.5 (X-X*)**2 (D2Y/DX2 )*  =  C(X**2) 

EX  =  -SQR  T ( ( Y (  1-1 l-YSTAR ) /C ) 


USE  NEWTON-RAPHSON  METHOD  TO  INVERT  QUARTIC  TO  OBTAIN  X  FROM  Y 
10  Kl=K0+2 
K2= I -1 

DO  50  K=K1,K2 
K3  =  K 1 +  K2-K 
WY=Y(K3  ) 

DO  20  L  =  1 » 1 00 

F  = (  (A*EX+B)*EX+C)*EX*EX+E-WY 

F 1 = ( <4.*A*EX+3.*B)*EX+2.*C)*EX 

ERR0R=-F/F1 

Q=ABS(ERROR)/EX 

EX=EX+ERROR 

IF  (Q. LT.. 000001 )  GO  TO  40 
20  CONTINUE 

WRITE  (6,30)  WY, EX, ERROR, Q 

30  FORMAT  (36H0  INLET  ITERATION  DID  NOT  CONVERGE  /1P4E15.7//) 
40  X ( K3  )  =  E X 
50  CONTINUE 
X ( I  )  =0. 

RETURN 

END 
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SUBROUTINF  SKIN 
C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT,XT,HO* 

1  DXDX*DRDX*DEMDX»DRH0DX*DPEDX*DUEDX*DHEDX, 

2  CF  *COEFF  »DELSTR*DELTA*DHIDX*DH1 *  DH2 ,DRWDX, 

3  'dsvdel,dsvth,dswdx,dthdx,emue,gas*had,hfi » 

4  HI *HJ*HK.»HW*IP*IP1,LAWCF*LM*M*N*PAGE,PR*Q» 

5  re»recvry»reth»reyanl*rhow»roe*rw»st*sumi , 

6  SUM2*TE»THETA,THVDEL,TW,TX(8)  ,  X  END , XPUNCH , 

7  XTRA1 *XTRA2*XTRA3*ZE 
DATA  AIR/6HAIR  / 

C 

c 

GO  TO  ( 10*60) ,  LAWCF 
C 
C 

C  SPAULDING  -  CHI  LAW  (REF.  J.FL.MECH.  JAN. 64) 

10  GO  TO  ( 20,30) ,  N 

20  FTHETA  = (HE/HW)**.702*( HAD/HW) **.772 
ASO=DH2/HW 

IF  (ASQ  .LT.  0.)  ASQ=0 • 

B=DH1/HW 

Dl  =  SORT <  4.*ASQ+B**2  ) 

D2  = ARS I N (  ( 2.*ASQ-B )/Dl)+ARSIN(B/Dl) 

FC=DH2/HE/D2**2 
IF  (CF.GT.O.)  GO  TO  30 

FCCF=2. /(5.77*ALOG10< FTHET A*R E TH ) +3 . 8 ) **2 
CF=FCCF/FC 

30  UGPLUS  =  SQRT (2./(FC*CF)  ) 

FTHRTH=FTHETA*RETH 
DO  40  INDEX=1»25 
X=0«4*UGPLUS 
XSQ=X**2 
EXPX=EXP ( X ) 

FUNC  =  -  FTHRTH  +  UGPLUS**2/6 . + ( ( 1 . -2 . /X ) *EXPX  +  2./X+1.  -XSQ*<30. 
1+X*(  1 5.+X*( 4.5+X )))/180.)/4.F 

DER I V  =  UGPLUS/3.+ (  ( ( XSQ  -2 . *X  +  2  .  ) *EXPX  -  2.1/XSQ  -X* ( 12 . +X* ( 9 • +X* 

1 { 3.6+X ) ) ) /36. ) / 1 2 . 

ERROR  =  -FUNC /DER I  V 
UGPLUS=UGPLUS+ERROR 

IF  ( ABS ( ERROR) /UGPLUS-. 0001 )  50*50»40 
40  CONTINUE 
50  FCCF=2. /UGPLUS**2 
CF=FCCF/FC 
RETURN 
C 
C 

C  LUDWIEG  -  TILLMANN  LAW  WITH  REFERENCE  ENTHALPY 

60  60  TO  (70*90)  »  N 
70  HREF=0.5*(HW+HE)+0.22*(HAD-HE) 

CALL  RHOTZ  ( GAS , PE , HREF , RHOREF  ,  TREF  , ZREF ) 

IF  (RHOREF. LE. 0.  )  RHOR EF= 3 . 5* P E/ ( 2 5 0 3 0 . *HR EF  ) 

IF  (GAS.EQ.AIR)  CALL  EMUAIR  (HREF, PE  *  EMUREF ) 

IF  (GAS. ME. AIR)  CALL  EMUN2  ( TREF  ,  RHOREF , EMUREF ) 

80  REREF=RHOREF*UE/EMUREF 

FACTOR=  0.246* ( RHOR EF/ROE ) ** 1.268 
90  CF=F ACTOR* 1 0 .**<-. 678*DSVTH ) / ( ( RE REF* THETA ) **0 . 268 ) 

RETURN 

FND 
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SUBROUTINE  SLOPE  < X^ , DELSTR ,DRDX , DRWDX ) 

REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 
RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

THIS  SUBROUTINE  COMPUTES  AN  APPROXIMATE  VALUE  OF  THE  WALL 
SLOPE  (DRW/DX).  THE  VALUE  OF  DRW/DX  IS  OBTAINED  BY  ADDING 
THE  CORE  SLOPE  (DR/DX)  TO  THE  X-DERIVATIVE  OF  THE  DISPLACEMENT 
THICKNESS  ( DELSTR ) •  D ( DELSTR ) /DX  IS  OBTAINED  BY  FITTING  A 
QUADRATIC  IN  X  TO  DELSTR  AT  THREE  CONSECUTIVE  POINTS. 

DATA  X2 »DX21 .DX32/-100. *0.  *0. / 

DRWDX=DRDX 
IF  ( X3-X2 )  10,30.20 
10  X2=-100. 

DX2 1  =  0  « 

DX32=0» 

20  OKAY  =  DX  21 
DX21=DX32 
DX32=X3-X2 
DR2 1 =DR  32 
X2  =  X3 
R2=R3 

DRDXP=DRDX1 
DRDX1=DRDX 
30  R3=DELSTR 
DR32=R3-R2 
IF  (OKAY)  40,50,40 

40  DRWDX2=(DR32*DX21/DX32+DR21*DX32/DX21)/(DX32+DX21 l+DRDXP 
DRWDX3=2.*DR32/DX32-< DRWDX2-DRDXP ) +DRDX 
DRWDX=DRWDX3 
50  RETURN 
END 


SUBROUTINE  TERM ( X , Z , D , F ) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

C  THIS  SUBROUTINE  PROVIDES  THE  TERMINATION  CONDITION  REQUIRED  BY 

C  FNOL2.  TERMINATION  WILL  BE  AT  X  =  XEND  UNLESS  XEND  EQUALS  OR 

C  EXCEEDS  999.  IN  WHICH  CASE  TERMINATION  WILL  BE  AT  NOZZLE  EXIT. 

C 

COMMON  EX,R,EM,RHOE,PE,UE,HE,REXIT,XEXIT,XT,HO, 

1  DXDX»DRDX,DEMDX»DRHODX,DPEDX»DUEDX,DHEDX, 

2  CF.COEFF, DELSTR, DELTA, DHIDX,DH1,DH2,DRWDX, 

3  DSVDEL,DSVTH»DSWDX,DTHDX,EMUE,GAS,HAD»HFI, 

4  HI ,HJ»HK,HW,IP,IP1,LAWCF»LM,M,N,PAGE,PR,Q, 

5  RE,RECVRY»RETH,REYANL,RHOW,ROE,RW,ST,SUMl , 

6  SUM2  »TE  , THETA, THVDEL  »TW  ,TX ( 8 )  , XEND , XPUNCH , 

7  XTRA1 »XTRA2»XTRA3»ZE 
DIMENSION  Z ( 5 )  ,D< 5  ) 

C 

IF  ( XEND. CE. 999. )  GO  TO  20 
C 

C  END  AT  SPECIFIED  LOCATION 

10  F=X-XEND 
GO  TO  30 
C 

C  END  AT  NOZZLE  EXIT 

20  F=X-XEXIT 
30  RETURN 
END 
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SUBROUTINE  WALL  (X) 

C 

C  REAL  GAS  TURBULENT  BOUNDARY  LAYER  PROGRAM 

C  RGTBLP  -  1/30/71  -  W.J.GLOWACKI 

C 

C 

C  THIS  SUBROUTINE  EVALUATES  THE  ENTHALPY  AND/OR  TEMPERATURE 

C  AT  THE  WALL  ACCORDING  TO  THE  OPTION  SPECIFIED  BY  M 

C 

C  IF  M=4 »  READ  TEMPERATURE  DISTRIBUTION  CARD 

C  IF  M=3 ,  COMPUTE  TW  FROM  ADIABATIC  WALL  ENTHALPY.  INITIAL 

C  GUESS  FOR  TW  AT  FIRST  POINT  MUST  BE  SUPPLIED. 

C  IF  M  =  2 ,  COMPUTE  TW  AND  HW  FOR  GIVEN  TEMPERATURE  DISTRIBUTION 

C  IF  M= 1 ,  TW  IS  CONSTANT  ALONG  NOZZLE  BUT  HW  VARIES  UNLESS  Z=1 

C 

COMMON  EX.R.EM.RHOE.PE.UE.HE.REXIT.XEXIT.XT.HO, 

1  DXDX .DRDX.DEMDX.DRHODX.DPEDX.DUEDX .DHEDX, 

2  CF.COEFF.DELSTR.DELTA.DHIDX.DHl ,DH2 .DRWDX, 

3  DSVDEL.DSVTH.DSWDX.DTHDX.EMUE.GAS.HAD.HFI . 

4  HI .HJ.HK.HW.IP.IPl.LAWCF.LM.M.N.PAGE.PR.Q. 

5  RE.RECVRY.RETH.REYANL.RHOW.ROE .RW.ST ,SUM1 » 

6  SUM2» TE, THETA, THVDEL , T W , T X < 8 )  , XEND  .XPUNCH , 

7  XTRA1.XTRA2.XTRA3.ZE 
C 

C 

GO  TO  ( 60.40.70,10)  ,  M 


READ  M,T*,X1»T1,X2,T2»TEXIT 
10  READ  (5,20  M ,  (  TX  (  I  )  ,  I  =2 , 6  )  ,  T  X  (  8  ) 
20  FORMAT  ( I5.6F10.5) 

TX ( 1 ) =0. 

TX(7)=XEXIT 
TW  =  TX ( 2  ) 

INITIAL  GUESS  FOR  RHOW 
RHOW=PE/ ( 3196. 8*TW) 

GO  TO  ( 80,30.80)  ,  M 


EVALUATE  CONSTANTS  OF  STATEMENT  30  FROM  INPUT  DATA 
30  DTOE=TX(2)-TX(8) 

F  =  AL0G(DT0E/(TX(4)-TX(8)  )  ) 

E  =  AL0G(DT0E/(TX(6)-TX(8)  )  ) 

D  =  TX ( 8  ) 

C  =  A  LOG 1 0 ( F/E) /ALOGlOt  TX( 3 ) /TX ( 5  )  ) 

B=F**(1./C)/TX(3) 

A=DTOE 
GO  TO  80 


EVALUATE  HW  FROM  TW 
40  IF  ( X.GE.O.  )  GO  TO  50 
TW= A+D 
GO  TO  60 

50  TW=A*EXP(-(B*X) **C ) +D 
60  CALL  HWALL  ( GAS , PE , TW , HW , RHOW  ) 
GO  TO  80 


EVALUATE  TW  FROM  ADIABATIC  HW 
70  HW=HAD 

CALL  RHOTZ  ( GAS , PE , HW , RHOW » TW , ZW ) 
80  RHOW=RHOW 
RETURN 
END 
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SUBROUTINE  FN0L2 ( J * N, G , L , M , NE ♦ X , V » D . OER I V , TERM .OUTPUT ) 
DIMENSION  Y  <  5  0 )  »  D ( 5  0  > ,  YB ( 30  » 6  )  »G I  2 ( 30 ) , G I  3 ( 30 ) »G 1 4( 30 ) *EF(30) 
1EF1 (30! *EF2(30) *EF3(30) *Y1(30) .ERROR (30) ,HA ( 30 ) . Y A ( 50 ) .DA ( 5 0 ) 
2 YC ( 30 )  *  YP ( 3  0  )  »  YD ! 50 ) 

DOUBLE  PRECISION  XD , YD , YA ♦ YC , YP , Y 1 
EC  =  Y ( N  +  3  ) 

1  H  =  G 

2  HZ  =H 

3  LN  =  N+M  A  XO  (  L  .  3  ) 

4  NA  =  0 

5  NB  =  1 

6  NF  =  0 

7  NG  =  0 

8  F  =  0. 

9  FA  =  0. 

10  FB=0. 

11  FC=0. 

12  FD=0 • 

13  ENE=NE 

DO  200  1  =  1  . LN 
200  YD (  I ) =DBLE ( Y  (  I  )  ) 

XD=DBLE ( X ) 

14  IF(J-3) 15,21.15 

15  IF(NE)18, 16,18 

16  J A  =  4 

17  GO  TO  22 

18  RE1=10.**(-ENE) 

19  RE2=10.**(-ENE-3.0) 

20  REM=10.**(-ENE-1.5) 

21  JA  =  1 

22  DO  25  1  =  1, N 

23  DO  24  I C  =  1 » 5 

24  YB ! I » IC ) =0. 

25  ERRORt I ) =0 • 

26  CALL  DERIV(X.Y.D) 

DO  300  1  =  1  , N 

G I  2  (  I  )  =  D  (  I  ) 

G  I  3 ( I ) =D ( I  ) 

G 1 4 ( I )=D( I ) 

300  EF( I )=D(  I) 

27  CALL  OUTPUT(X,Y,D,  ERROR,  N.L.H) 

28  FD=Y (  N  +  l  ) 

29  I F ( J-2 )  30,129,30 

30  GO  TO(31, 37, 35,37), JA 

31  DO  33  1=1, LN 

32  YA ( I ) = YD ( I ) 

33  DA  < I  >  =D (  I  ) 

34  60  TO  37 

35  HB=H 

36  H=2«*H 

37  HD2  =  .  5*H 
DO  39  1  =  1, N 

38  YB ( I ,NB)=D( I ) 

XL  =  D( I)  *  HD2 

39  Y ( I ) =SNGL ( YD ( I )+XL) 

X=SNGL ( XD+HD2 ) 

40  CALL  DERIV  (X,Y,GI2) 

41  DO  42  1=1, N 

XL  =  G I  2 (  I  ) *HD2 

42  Y ( I ) =SNGL ( YD (  I  )+XL) 

43  CALL  DERIV  (X,Y,GI3) 

44  DO  45  1  =  1, N 
XL  =  G  I  3 (  I )*H 

45  Y (  I  ) =SNGL ( YD ( I )+XL) 

X=SNGL ( XD+H ) 

46  CALL  DERIV(X,Y,GI4) 
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47  HD6  =H/6.  X 

GO  TO(48*55*60»66)*JA 

48  DO  52  I = 1  * N 

XL  = ( D ( I  )  +  2 • * ( G I  2 (  I  )  +  G I  3 (  I  )  )  +GI4(1))*HD6 

49  YC ( I ) =  YD ( I ) +XL 

51  YD ( I ) =Y A (  I  ) 

52  ERROR! I ) =0 . 

53  JA=3 

54  GO  TO  35 

55  DO  57  1=1, N 

XL=(D(I)  +  2 • * ( G I  2 ( I )  +  G  I  3 ( I  )  )  +GI4(I))*HD6 

56  YD! I ) = YD ( I ) +XL 

57  ERROR! I ) =SNGL( YD!  I  ) -YP ( I )  ) / 1 5  • 

58  JA= 1 

59  GO  TO  681 

60  DO  62  1=1, N 

61  YD! 1 ) = YC ( I  ) 

X  L  =  (  D  (  I  )  +  2  •  *  (  G  I  2  (  I  )  +'  G  I  3  (  I  )  )  +GI4(I))*HD6 

62  YP! I ) = Y A ( I  ) +XL 

63  H=HB 

64  JA  =  2 

65  GO  TO  681 

66  DO  68  1  =  1, N 

XL= ( D ( I  )  +  2.*(GI2(I)  +  G  I  3 ( I  )  )  +GI4(I))*HD6 

67  YD! I ) = YD ( I )+XL 

68  ERROR ( I ) =0. 

681  DO  69  1  =  1, N 

69  Y (  I ) =  SNGL ( YD ( I )  ) 

XD=XD+H 

X=SNGL ( XD ) 

70  CALL  DE R I V ( X , Y , D ) 

71  FC=  F 

72  CALL  TERM(X,Y,D,F) 

73  IF ( ABS! F1-1.0E-5  1731,731,733 

731  NF=5 

732  GO  TO  124 

733  IF(F)74»124,76 

74  FA  =  1  • 

75  GO  TO  77 

76  FB=  1  • 

77  IF(FA-FB)83,78,83 

78  NF=NF+1 

79  JA=4 

80  NB  =  1 

81  H=H*F/ ( FC-F 1 

82  IF(NF-4)37, 37*124 

83  IF! NE 184, 117,84 

84  IF! JA-1 1117,85, 117 

85  I F ( J-3 186,117,86 

86  DO  95  1=1, N 

I F ( Y ! I  1  1886,885,886 

885  HA! I  1=1000. 

GO  TO  95 

886  I  F ( EC  1 880 , 890 , 87 

87  I F ( ABS ( Y ( I  1  1 -EC  1  880,880,890 

880  IF(ABS(ERROR( I  1 1  —  R  E  2 1  882,94,881 

881  IF (ABS! ERROR! I  1  1 -RE  1 1 94 , 9 4 , 88 2 

882  HA!  I  1 =H*( REM/ (ABS (ERROR! I  1 1 +. 0000000001 1 >**< .2  1 

883  GO  TO  95 

890  IF! ABS! ERROR!  I  1 / Y ( I  1 1 -R E2 1 8 92 , 94 , 89 1 

891  IF ( ABS ( ERROR ( 1  1 / Y ( I  1 1 -R E 1 1 94 , 94 , 8 92 

892  HA!  I  )=H*( REM/ (ABS (ERROR!  I  1 / Y (  I  1  1 +  . 0000000001 1  1 ** (  .  2 1 

893  GO  TO  95 

94  HA ( I  1 =H 

95  CONTINUE 

96  HB=HA ( N 1 
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97  66  98  I=1»N 

98  HB=AMI N 1 ( HA ( I ) *HB ) 

99  IF(H— HB)100*117*10l 

100  IF(HZ-H)101*101*116 

101  DO  103  I =1 *LN 

102  YD ( I ) =Y A ( I  ) 

Y ( I ) =SNGL ( YD ( I ) ) 

103  D  < I > =DA  (  I  ) 

104  IF1NB-6 1107*105*105 

105  XD=XD-H 

106  GO  TO  109 

107  XD=XD-2.*H 

108  HZ  =H 

109  H=HB 
X=SNGL ( XD) 

CALL  DERIV(X*Y*D) 

110  NB=1 

111  XABS=ABS( ,000001*X ) 

112  IF(ABS(H)-XABS)113*113*117 

113  NG=NG+ 1 

114  H=SIGN( XABS.HB) 

115  I F ( NG-10)  124*126  *  126 

116  HZ  =H 

117  IF(M) 118*118*121 

118  I F ( ABS  (  Y(N  +  l)-FD)-Y(N+2 )) 29 ,1 19, 119 

119  FD= Y (  N  +  l  ) 

120  GO  TO  124 

121  NA=NA+1 

122  IF(M-NA) 123,123,29 

123  NA  =  0 

124  CALL  OUTPUT (X,Y,D  *  ERROR , N  »  L  ,H  ) 

125  IF(NF-4)29»29*126 

126  WRITE  (6,127) 

127  FORMAT ( 1H0 ) 

128  RETURN 

129  NB=NB+ 1 

130  IF(NB-6)30»131,136 

131  DO  134  1=1* N 

132  EF3 ( I >=YB< 1,3) 

133  EF2 ( I )  =  YB  ( I ,4) 

134  EF1 ( I )  =  YB ( 1,5) 

135  GO  TO  137 

136  NB= 1 0 

137  HD24  =H/24. 

DO  138  1=1,  N 

XL  =(55.*D(I)  — 59 • *  EF 1 (  I  )  +37.*EF2(I)  -9.*EF3 ( I ) ) *HD24 
YP (  I ) = Y  D ( I ) +XL 

138  Y ( I ) =SNGL ( YP (  I  )  ) 

X  =  SNGL ( XD+H ) 

139  CALL  DERIV(X,Y,EF) 

140  DO  142  1=1, LN 

141  YA( I ) = YD ( I ) 

142  DA ( I ) =D ( I ) 

143  DO  148  1  =  1, N 

XL  =  <9.*EF(I)  +19.*D(I)  -5.*EF1U)  +EF2(I))*HD24 

144  YD (  I ) = Y D ( I ) +XL 

145  ERROR ( I ) =-SNGL ( YD ( I ) -YP ( I ) ) / 1 4 . 

146  EF3( I )  =  E F 2 (  I ) 

147  EF2 (  I ) =  EF 1 (  I  ) 

148  EF1 (I )=D( I ) 

149  GO  TO  681 

END 
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